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1  INTRODUCTION 


Nascap-2k  is  a  spacecraft  charging  and  plasma  interactions  code  designed  to  be  used  by 
spacecraft  designers,  aerospace  and  materials  engineers,  and  space  plasma  environments  experts 
to  study  the  effects  of  both  the  natural  and  spacecraft-generated  plasma  environments  on 
spacecraft  systems. 

Designers  of  spacecraft  for  government,  commercial,  and  research  purposes  require  advanced 
modeling  capabilities  to  guide  the  design  of  satellites  that  can  survive  and  operate  properly  in  the 
natural  environment.  Computer  modeling  of  flight  experiments  (such  as  SCATHA  (Spacecraft 
Charging  at  High  Altitude),  the  SPEAR  (Space  Power  Experiments  Aboard  Rockets)-’  series, 
and  CHAWS  (Charging  Hazards  and  Wake  Studies)4)  has  demonstrated  excellent  ability  to 
predict  both  steady-state  and  dynamic  interactions  between  high-voltage  spacecraft  and  the 
ambient  plasma.  This  ability  was  also  extended  to  inherently  dynamic  problems  involving  three- 
dimensional  space  charge  sheath  formation,  current  flow  in  the  quasi-neutral  presheath, 
breakdown  phenomena,  plasma  kinetics,  ionization  processes,  and  the  effect  of  unsteady 
processes  on  spacecraft  charging. 

Nascap-2k  encapsulates  the  knowledge  gained  in  these  efforts  as  well  as  in  modeling  spacecraft 
with  unique  requirements  and  destined  for  disparate  environments.  This  gives  the  spacecraft 
designer  quality  modeling  capabilities  by  taking  advantage  of  the  present  understanding  of  the 
pertinent  phenomena,  employing  advanced  algorithms,  and  implementing  a  state-of-the-art  user 
interface,  including  three-dimensional  post-processing  graphics. 

The  core  capabilities  of  Nascap-2k  are  as  follows: 

■  Define  spacecraft  surfaces  and  geometry  and  the  structure  of  the  computational  space 
surrounding  the  spacecraft. 

■  Solve  for  time-dependent  potentials  on  spacecraft  surfaces. 

■  Solve  the  electrostatic  potential  around  the  object,  with  flexible  boundary  conditions  on  the 
object  and  with  space-charge  computed  either  fully  by  particles,  fully  analytically,  or  in  a 
hybrid  manner. 

■  Generate,  track,  and  otherwise  process  particles  of  various  species,  represented  as  macro¬ 
particles,  in  the  computational  space. 

■  View  surface  potentials,  space  potentials,  particle  trajectories,  and  time-dependent  potentials 
and  currents. 

Nascap-2k  calculates  surface  charging  in  tenuous  plasma  environments,  such  as  geosynchronous 
earth  orbit  (GEO)  and  interplanetary  (Solar  Wind)  environments,  using  the  Boundary  Element 
Method5  (BEM).  The  Boundary  Element  Method  facilitates  calculation  of  surface  electric  fields 
(which  limit  the  emission  of  photoelectrons  and  secondary  electrons)  without  the  need  to  grid  the 
space  surrounding  the  spacecraft.  It  also  enables  Nascap-2k  to  anticipate  electric  field  changes 
due  to  surface  charging,  resulting  in  a  smoother  and  more  stable  charging  simulation. 

Nascap-2k  uses  a  high-order,  finite-element  representation  for  the  electrostatic  potential  that 
ensures  electric  fields  are  strictly  continuous  throughout  space.  The  electrostatic  potential  solver 
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uses  a  finite  element/conjugate  gradient  technique  to  solve  for  the  potentials  and  fields  on  the 
spacecraft  surface  and  throughout  the  surrounding  space.  Space  charge  density  models  presently 
include  Laplacian,  Linear,  Nonlinear,  Frozen  Ions,  Consistent  with  Ion  Density,  Full  PIC 
(Particle-in-Cell),  and  Hybrid  PIC. 

Particle  tracking  is  used  to  simulate  sheath  currents,  to  study  detector  response,  or  to  generate 
space  charge  evolution  for  dynamic  calculations.  Nascap-2k  generates  macroparticles  (each  of 
which  represents  a  collection  of  particles)  at  a  sheath  boundary,  the  problem  boundary,  specified 
surface  elements,  throughout  all  space,  or  at  user-specified  locations.  Particles  are  tracked  for  a 
specified  amount  of  time,  with  the  timestep  automatically  subdivided  at  each  step  of  each 
particle  to  maintain  accuracy.  The  current  to  each  surface  element  of  the  spacecraft  and 
optionally  through  each  volume  element  is  recorded  for  further  processing. 

Nascap-2k  User  s  Manual  is  the  primary  document  on  the  use  of  the  code.  This  document, 
Nascap-2k  Scientific  Documentation,  describes  the  physics  and  numeric  models  used  in  the 
surface  charging,  potential  solution,  and  particle  tracking  portions  of  the  code.  Section  2, 
Physical  Models,  describes  the  physics  models  included  in  the  code.  Section  3,  Numeric  Models, 
describes  the  implementation. 

2  PHYSICAL  MODELS 

Nascap-2k  incorporates  physical  models  appropriate  to  both  tenuous  (e.g.,  GEO  orbit  or 
interplanetary  missions)  and  dense  (e.g.,  LEO  orbit)  plasma  environments.  The  code  solves  for 
environmentally-induced  time-dependent  potentials  on  spacecraft  surfaces,  the  electrostatic 
potential  about  the  object,  and  the  resulting  charged  particle  motion. 

Time-dependent  surface  currents  are  computed  using  orbit  limited  currents  in  tenuous  plasma, 
tracked  currents  and/or  analytic  approximations  in  dense  plasma,  or  a  combination  in 
intennediate  environments. 

The  electrostatic  potential  solver  uses  a  conjugate  gradient  technique  with  flexible  boundary 
conditions  to  solve  Poisson’s  equation  for  the  potentials  and  fields  throughout  the  computational 
space.  Several  analytic  and  numeric  space  charge  density  models  are  available. 

Particle  tracking  is  used  to  study  sheath  currents,  to  study  detector  response,  to  generate  steady- 
state  charge  densities,  and  to  generate  space  charge  evolution  for  dynamic  calculations. 

2,1  Surface  Charging  from  Orbit  Limited  Currents 

Spacecraft  surface  charging  is  the  accumulation  of  charge  on  spacecraft  surfaces.  As  illustrated 
in  Figure  1,  several  different  current  components  can  contribute  to  the  charging.  ’  ’  ’  First  there 
are  the  incident  charged  particles.  The  impact  of  high-energy  electrons  and  ions  causes  the 
ejection  of  secondary  electrons  into  space.  Incident  electrons  can  also  be  reflected  as  backscatter. 
In  sunlight,  solar  ultraviolet  radiation  generates  low-energy  photoelectrons  which  are  ejected 
from  the  surface.  Due  to  the  low  mass  of  electrons,  the  incident  electron  current  is  generally 
higher  than  the  ion  current.  By  itself,  this  would  lead  to  negative  charging.  However,  with  the 
inclusion  of  the  other  terms,  all  of  which  contribute  positive  current,  the  net  charge  can  be  either 
positive  or  negative. 
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Secondary,  Backscattered,  Photo  Electrons 

Figure  1.  High  Negative  Potentials  Can  Result  from  the  Accumulation  of  Charge  on  Spacecraft 
Surfaces. 

The  rest  of  the  spacecraft  influences  the  potential  and  local  electric  field  of  each  surface.  Each 
conductor  and  each  insulating  spacecraft  surface  collects  current  from  the  plasma  and  is 
capacitively  and  resistively  coupled  to  the  other  conductors  and  other  surfaces.  Additionally, 
some  spacecraft  have  current  sources,  such  as  electron  and  ion  guns,  that  contribute  to  the  net 
current.  Nascap-2k  uses  the  spacecraft  geometry,  surface  materials,  plasma  environment,  and  sun 
intensity  and  direction  to  calculate  the  time  history  of  the  surface  potentials,  electric  fields,  and 
fluxes. 

If  the  net  current  is  initially  negative  or  positive,  the  exposed  surface  begins  to  acquire  a  negative 
or  positive  potential  respectively.  As  the  magnitude  of  the  potential  increases,  the  net  current  is 
attenuated,  until  it  eventually  approaches  zero  and  the  surface  potential  reaches  its  equilibrium 
value.  Negative  potentials  as  high  as  ten  kilovolts  have  been  observed  during  geosynchronous 
substorms  and  over  one  kilovolt  during  auroral  precipitation  events.  Because  positive  currents 
generally  result  from  emission  of  low  energy  secondary  electrons  and  photoelectrons,  the 
positive  potentials  that  can  be  attained  are  relatively  modest.  As  soon  as  the  surface  reaches  a 
potential  greater  than  that  of  the  emitted  electrons  they  are  re-attracted  to  the  surface,  so  their 
contribution  to  the  current  is  suppressed.  The  equilibrium  potential  is  determined  by  the 
suppression  of  low  energy  electron  emission  due  to  the  surface’s  own  electric  field.  A  similar 
suppression  effect  may  occur  due  to  the  electric  fields  of  neighboring  negative  charged  surfaces. 
This  effect  is  the  primary  reason  that  accurate  spacecraft  surface  charging  calculations  are 
inherently  three  dimensional. 

Positive  potentials  up  to  100  V  are  possible  in  environments  where  the  plasma  current  is  small 
compared  with  the  photo  current  (outer  magnetosphere  and  in  the  solar  wind).  In  these 
environments,  the  small  incident  plasma  current  is  balanced  by  escape  of  the  high  energy  tail  of 
the  photoelectron  energy  distribution,  even  though  the  bulk  of  the  emitted  photoelectrons  return 
to  the  surface. 

Nascap-2k  models  the  spacecraft  as  a  collection  of  surfaces  that  are  either  conductive  or  are  a 
dielectric  film  over  a  conductive  substrate.  The  conductive  surfaces  and  substrates  can  be  at  the 
same  potential  or  biased  with  respect  to  each  other. 
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2.1.1  Plasma  Environment 

In  the  orbit  limited  approximation  of  the  current  incident  on  a  surface,  the  flux  is  calculated  as 
that  incident  on  a  sphere  of  the  same  potential  and  depends  only  on  the  surface  potential,  and 
optionally  an  angle  with  respect  to  the  direction  of  flow.  This  approximation  applies  if  any 
charged  particle  far  from  the  spacecraft  can  reach  the  surface,  that  is,  there  are  no  excluded 
orbits.  For  spacecraft  this  is  generally  true  if  the  Debye  length  is  long  compared  with  the 
spacecraft  size  and  the  spacecraft  is  essentially  convex  with  no  self-shading.  Calculations  using 
the  “Surface  Charging”  with  “Analytic  Currents”  option  use  this  approximation  for  both  species, 
and  calculations  using  the  “Surface  Charging”  with  “Tracked  Ions  and  Analytic  Electron 
Currents”  option  use  this  approximation  for  electron  currents. 

The  analytic  current  density  to  a  surface  at  potential  (j>  is  given  by  integrating  over  the 
distribution  function. 


oo  oo 


j  =  q 

2Ee2 
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Where  the  integration  variable  E  represents  the  energy  at  the  surface,  the  integration  limit  L  is  0 
for  the  repelled  species  and  |(j)|  for  the  attracted  species,  the  fraction  in  parentheses  represents  the 
orbit  limited  enhancement  or  reduction  of  current,  the  upper  sign  is  for  ions  and  the  lower  sign  is 
for  electrons,  f(E)  is  the  distribution  function  at  infinity,  and  F(E)  is  the  flux  at  infinity. 
Depending  on  the  model  selected,  the  flux  is  given  by  one  of  the  following  formulas  in  which  the 
particle  energy  E  and  the  temperature  9  are  measured  in  electron  volts,  and  the  surface  potential 
cj>  is  measured  in  volts. 
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This  distribution  is  used  to  model  auroral  electrons. 
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where  H  is  the  Heaviside  step  function  and  El  and  Eu  are  upper  and  lower  energy  limits. 

Kappa 


F(E) 


27lK0m  K  0 


-n 


r(*+i) 

rfK- 


Convected  Maxwellian 

This  distribution  is  used  to  model  a  flowing  plasma. 

E  +  mU'  /  2  -  sj  Em/  2U  cos  x 
0 

v  J 


F(E,x)  = 


2710m  0 


nexp 


(5) 


(6) 


where  U  is  the  velocity  of  the  plasma  in  the  spacecraft  frame  of  reference  and  %  is  the  angle 
between  the  flow  vector  and  the  incident  velocity  at  infinity,  which  can  be  related  to  the  velocity 
at  which  the  particle  strikes  the  surface.  In  this  case,  obtaining  the  flux  to  the  surface  requires 
integrating  over  angle  as  well  as  energy. 

Measured 

This  option  should  be  used  with  extreme  caution  as  inadequately  characterized  environments  can 
lead  to  numeric  instabilities  or  unphysical  results. 


F(E)  =  £  fj  (E)  where  f;  (E) 
i 


Fj  forEj  <  E  <  Ei+1 
0  otherwise 


(V) 


2,1.2  Secondary  and  Backscattered  Current 

The  secondary  and  backscattered  currents  are  given  by 


Js,b 


Ys’b(E) 


E  ± 


F(E  ±  t())dE 


(8) 


where  Ys,b  is  the  number  of  secondary  or  backscattered  (respectively)  electrons  per  incident 
charged  particle  (after  accounting  for  the  angular  distribution  of  the  incident  particles)  and  F  is 
the  flux  at  infinity  of  the  appropriate  incident  charged  particle. 


We  take  the  spectrum  of  emitted  secondary  electrons  to  be  a  Maxwellian  characterized  by  a 
temperature  of  2  eV.  It  follows  that  the  effective  secondary  electron  current  from  surfaces  at 
positive  potentials  is  reduced  by  a  factor  of  exp(-(j)/2)  relative  to  Equation  (8).  If  the  surface 
potential  is  negative,  the  secondary  electron  current  can  still  be  reduced  by  an  electron-attracting 
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(E  n  >  0)  electric  field  caused  by  other,  more  negative  surfaces,  by  negative  potentials  in  a 
nearby  wake,  or  by  the  space  charge  of  the  emitted  electrons  themselves. 

2.1.3  Photocurrent 

For  surfaces  at  negative  potentials,  the  photocurrent  depends  on  the  surface  material,  angle  of 
incident  sunlight,  and  distance  from  the  sun.  As  with  secondaries,  the  photoelectron  spectrum  is 
by  default  taken  to  be  a  Maxwellian  characterized  by  a  temperature  of  2  eV,  so  that  current 
escaping  from  surfaces  at  positive  potentials  is  reduced  by  a  factor  of  exp(-(j)/2). 

Spacecraft  in  the  solar  wind  normally  charge  to  positive  potential  (a  few  volts  to  tens  of  volts)  to 
reduce  the  emitted  photoelectron  current  to  balance  a  very  low  incident  current  of  ambient 
electrons.  Unlike  geosynchronous  substonn  charging,  where  the  spacecraft  goes  negative, 
detennining  how  positive  the  spacecraft  charges  requires  knowledge  of  the  high-energy  portion 
of  the  photoelectron  spectrum.  To  override  the  default,  the  user  can  supply  this  spectrum  as  either 
a  table  or  as  the  sum  of  exponentials  of  the  form  Ae  1  0 . 

In  magnetospheric  or  interplanetary  space  the  space  charge  near  a  photoemitting  surface  is 
dominated  by  photoelectrons.  This  causes  the  potential  to  decrease  with  distance  from  the 
surface  until  loss  of  returning  electrons  together  with  spatial  divergence  reduces  the 
photoelectron  density  to  the  order  of  the  ambient  plasma  density.  Both  close  (~0. 1  AU)  to  the  sun 
where  photoemission  is  high  and  on  a  very  large  spacecraft  (e.g.,  Solar  Sail)  for  which 
divergence  is  low,  the  potential  some  distance  from  the  surface  can  reach  negative  values.  When 
this  is  the  case,  the  minimum  energy  for  photoelectron  escape  is  increased  beyond  the  surface 
potential.  The  barrier  reduces  the  net  photocurrent  and  lowers  the  surface  potential.  This  problem 
was  studied  analytically  by  Guernsey  and  Fu, 10  and,  more  recently,  numerically  by  Ergun  et  al. 1 1 

When  the  user  requests  that  this  space  charge  limiting  of  photocurrent  be  included  in  the  surface 
potential  calculation,  the  code  computes  a  barrier  height.  If  the  space  charge  barrier  exceeds  the 
surface  potential,  the  code  reduces  the  net  photocurrent  and  the  secondary  electron  current.  In 
addition  to  the  photoemission  current  and  the  ambient  plasma  density,  the  preliminary  Nascap-2k 
model  depends  on  the  effective  radius  of  curvature  of  the  surface  element  (estimated  by  (j)/E, 
where  E  is  the  electric  field  [calculated  using  the  boundary  element  method  for  the  spacecraft  at 
uniform  potential  cj)]  and  constrained  to  be  between  0.3  m  to  2.0  m)  and  the  shape  of  the 

photoemission  spectrum  (taken  to  be  J  =  53e~E/1'6  +  21e~E/3  +4e^E/8'9  12).  The  effect  of  this 
algorithm  on  large,  flat  surfaces  is  to  reduce  the  effective  photoemission  current  from  the  interior 
of  the  surface  (where  effective  radius  of  curvature  is  large)  but  not  near  the  edges  (where  the 
effective  radius  of  curvature  is  small). 

2.1.4  Other  Current  Sources 
Conduction  Current  through  Insulators 

In  surface  charging  calculations,  the  charge  deposition  thickness  is  short  compared  with  the 
thickness  of  the  insulator  layer,  and  therefore,  the  electric  field  is  spatially  uniform  throughout 
the  layer  and  current  is  conducted  through  it.  (Note:  the  range  of  a  40  keV  electron  in  aluminum 
is  ~15  microns  [~  0.0006  inches].)  For  insulating  surface  elements,  the  current  through  the  layer 
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to  the  underlying  conductive  substrate  is  given  as  a  function  of  the  conductivity  and  thickness  of 
the  layer  by  Ohm’s  Law: 


J  conduction 


A(J) 

pd 


(9) 


where  A<j)  is  the  difference  in  potential  between  the  surface  and  the  conductive  substrate  and  p 
and  d  are  the  resistivity  and  thickness  of  the  layer  respectively. 

Surface  Conductivity 

Nas cap-2 k  accounts  for  current  flowing  along  surfaces.  Current  flows  between  insulating  surface 
elements  of  a  common  material  and  with  a  common  edge,  transporting  charge  over  a  wide 
expanse  of  such  material.  The  surface  element  to  surface  element  conductivity  is  proportional  to 
the  length  of  the  common  edge  and  inversely  proportional  to  the  sum  of  the  centroid-to-edges 
distances.  For  a  surface  element  i  having  a  common  edge  with  j,  R,  =  (d/l,,)  p,  where  p  is  the 
surface  resistivity,  dj  is  the  distance  from  the  centroid  of  surface  element  i  to  the  center  of  the 
common  edge,  and  l,,  is  the  length  of  the  common  edge.  The  surface  conductance  across  the  edge 
between  adjoining  surface  elements  i  and  j  is 


1 

C  onductivity  =  — - — —  .  (10) 

Ri+Rj 

A  surface  can  be  grounded  by  a  strip  at  an  element  edge  (grounding  edge)  or  by  a  circular  contact 
located  at  a  node  (grounding  node). 

Charged  Particle  Emission 

Nascap-2k  also  allows  for  the  inclusion  of  the  current  due  to  charged  particle  emission  from  a 
conductor  (e.g.,  an  electron  gun).  This  capability  is  accessed  by  editing  the  script. 

2,1.5  Timescales 


Charged  particles  with  energies  below  50  keV  penetrate  less  than  a  hundred  microns  into  the 
spacecraft  skin.  While  the  time  for  the  entire  spacecraft  to  achieve  net  current  balance  is  typically 
several  milliseconds,  the  time  for  each  surface  to  achieve  its  own  equilibrium  potential  is 
thousands  of  times  longer.  The  development  of  differences  between  the  potentials  of  different 
surfaces  is  referred  to  as  differential  charging. 

The  rate  of  absolute  charging  is  determined  by  the  capacitance  of  the  spacecraft  to  infinity,  while 
the  rate  of  differential  charging  is  determined  by  the  capacitance  across  the  insulator  layer. 
Nascap-2k  does  not  take  the  presence  of  a  plasma  into  account  when  computing  the  capacitance 
to  space.  The  rate  of  change  of  the  potential  on  a  sphere  is  given  by 


d(|)  _  1  I  _  JR 

dt  4tts0  R  s0 


(11) 
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where  I  is  the  current  to  the  sphere,  J  is  the  current  density  to  the  sphere  surface,  and  R  is  the 
sphere  radius.  For  a  1 -meter  sphere  and  a  net  current  of  1  pA  m'  ,  the  sphere  begins  to  charge  at  a 
rate  of  100  kV  s'1.  As  the  sphere  charges  toward  its  steady-state  (zero  net  current)  potential,  the 
net  current  decreases  toward  zero  as  electrons  are  repelled.  Therefore,  the  charging  rate  decreases 
with  time. 


The  rate  of  change  of  the  difference  in  potential  across  a  dielectric  layer  is  given  by 

d(j)  _  Jd 
dt  s0 


(12) 


where  d  is  the  thickness  of  the  layer.  A  100  pm  layer  with  the  same  incident  current  density 
charges  at  a  rate  of  10  V  s'1.  In  this  representative  example,  differential  charging  takes  place  104 
times  more  slowly  than  absolute  charging. 

2.1.6  Shadowing 

A  surface  element  is  shadowed  from  the  sun  if,  either 

(1)  Its  nonnal  points  anti-sunward,  or 

(2)  Its  centroid  lies  behind  the  projection  of  another  sun- facing  surface. 

Shadowing  from  the  ram  flow  is  handled  in  the  same  manner. 

2.1.7  Circuit  Model 


Figure  2  shows  a  circuit  diagram  for  a  spacecraft  with  one  insulating  surface  element  and 
exposed  conducting  surfaces.  The  widely  differing  capacitances  of  surfaces  to  infinity,  Ca  and  Cs 
and  of  the  surface  to  spacecraft  ground,  CAs,  make  this  a  complex  numeric  problem. 

CAS  =8£0-~SxlO-6  Farad  (13) 

d 


C, 


Cs  «  47ts0  R 


4ttR2 


r  c  x 


vRy 


xlO  Farad 


(14) 


where  s,  d,  and  S  are  the  dielectric  constant,  thickness  (m),  and  surface  area  (m  )  of  the 
insulating  surface  element,  and  R  is  the  radius  of  a  sphere  approximating  the  size  of  the 
spacecraft.  The  potentials  as  a  function  of  time  are  computed  using  implicit  time  integration  of 
the  charging  equations,  which  relate  the  derivative  of  the  potential,  ®,  with  time  to  the  current,  I. 

cA®A+cAS(®A-®s)=iA 

-cas(®a-®s)+cs®a=ia 
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Insulating  surface 


All  conductive  surfaces 


Figure  2.  Circuit  Model  of  a  Spacecraft  with  One  Insulating  Surface  Element. 

The  multi-surface  problem  is  solved  by  linearizing  the  currents  and  inverting  the  matrix. 

C0>  =  l(0>)  (15) 

2.2  Material  Properties 

The  properties  of  the  spacecraft  surface  materials  determine  the  charging  rates  and  equilibrium 
surface  potentials.  Nascap-2ku  ses  the  models  of  secondary  electrons  due  to  incident  electrons, 
secondary  electrons  due  to  incident  ions,  backscattered  electrons,  and  photoelectrons  originally 
developed  for  NASCAP/GEO. 13 

2.2.1  Secondary  Electron  Emission  Due  to  Electron  Impact 

Secondary  electrons  are  those  emitted  from  a  surface  with  energies  below  50  eV  in  response  to 
energy  deposited  near  the  surface  by  incident  electrons.  Their  energy  distribution  is  usually 
peaked  below  10  eV.  The  secondary  yield,  Yee,  is  the  ratio  of  primary  to  secondary  electron 
current. 


emitted  secondary  current  due  to  electron  impact 
primary  electron  current 


(16) 


A  typical  curve  is  shown  in  Figure  3. 
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Figure  3.  Electron  Secondary  Yield  as  a  Function  of  Incident  Energy. 


The  secondary  electron  emission  yield  can  be  calculated  using  the  empirical  formula14: 


Yf 


ee 


W  =  c}( 


R 


dE 


0  dx 


-ax  cos 


vdx 


(17) 


where  x  is  the  path  length  of  penetration  of  a  primary  electron  beam  into  the  material,  R  is  the 
“range,”  or  maximum  penetration  length,  and  \j /  is  the  angle  of  incidence  of  the  primary  electron. 

This  equation  is  based  on  a  simple  physical  model 15 : 

1 .  The  number  of  secondary  electrons  produced  by  the  primary  beam  at  a  distance  x  is 
proportional  to  the  energy  loss  of  the  beam  or  “stopping  power”  of  the  material, 

dE 
dx  ' 

2.  The  fraction  of  the  secondaries  that  migrate  to  the  surface  and  escape  decreases 

exponentially  with  depth  (f  =  e_ax  cos  V|/  ).  Thus  only  secondaries  produced  within  a  few 
multiples  of  the  distance  1/a  (the  depth  of  escape)  from  the  surface  contribute  significantly 
to  the  observed  yield. 


The  stopping  power  for  incident  electrons  of  energy  E  is  related  to  the  range  of  these  electrons 
through  the  equation 


S(E) 


dE 

dR 

dx 

dE 

(18) 


The  usual  fonnulation  for  the  range  is  that  it  increases  with  the  energy,  E  of  the  incident 
electrons  in  a  way  that  approximates  a  simple  power  law16: 

R  =  bEq  (19) 
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where  q  is  a  bit  less  than  2  for  electrons  with  energy  <  300  keV.  Because  the  primary  beam  loses 
energy  as  it  passes  through  the  material,  E  and  S(E0,  x)  depend  on  the  path  length  x,  where  E0  is 
the  initial  electron  energy.  The  stopping  power  can  be  written  as 


S(E0,x) 


dE 

= 

dR 

1 

r  b  i 

dx 

dE 

qb 

Lr-xJ 

(20) 


Figure  4  shows  S(E0,  x)  plotted  against  x  for  several  values  of  E0.  Inspection  of  Figure  4  and  the 
equation  for  S(x)  illustrates  the  following  points: 

1 .  S(E0,  x)  increases  with  x,  slowly  at  first,  before  reaching  a  singularity  as  x  approaches  R. 

2.  The  initial  value  of  S(E0,  x)  decreases  with  increasing  initial  energy  E0. 

Both  of  these  observations  are  due  to  the  decrease  in  electron-atom  collision  cross-section  with 
increasing  energy. 


Path  Length,  x 


* —  1/a  — ► 
Depth  of  Escape 


(a) 


(b) 


Figure  4.  Energy  Deposition  Profiles  of  Normally  Incident  Primary  Electrons  for  Four  Incident 
Energies  E°  ,  ,  E^  ,  and  E^  and  Corresponding  Yield  Curve. 

The  yield  is  only  sensitive  to  the  details  of  the  stopping-power  depth-dependence  for  initial 
energies  with  ranges  of  the  same  order  as  the  escape  depth,  R  ~  1/a  (i.e.,  about  the  maximum  of 
the  yield  curve).  For  lower  energies,  R  «  1/a,  essentially  all  the  primary  energy  is  available  for 
detectable  secondary  production,  leading  to  a  linear  increase  in  yield  with  increasing  E0.  At 
higher  energies,  where  R  »  1/a,  S(E0,x)  remains  almost  constant  over  the  depth  of  escape. 
Therefore,  along  with  S(E0,x),  the  yield  decreases  as  E0  increases. 

Taking  this  into  account,  the  stopping  power  can  be  approximated  by  a  linear  expansion  in  x, 
about  x  =  0. 
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dE 

dx 


x 


(21) 
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The  simple  power  law  does  not  adequately  describe  the  available  experimental  data.  The  range 
has  different  exponents  for  low  energy  and  high  energy.  A  bi-exponential  range  law  with  four 
parameters  bi,  b2,  qi,  q2  fits  the  experimental  data  better. 

R=b!E^  +  b2Eg2  (22) 


Figure  5  shows  CSDA  (Continuous  Slowing  Down  Approximation)  range  values  for  Aluminum 
for  energies  from  20  eV  to  100  keV.  Note  the  drastic  change  in  behavior  near  300  eV.  The  bi¬ 
exponential  formula  R  =  1 42. 63E 0,2335  +  245. 15E1  7269  fits  well  throughout  the  range  plotted. 


Electron  Energy  (keV) 


Figure  5.  CSDA  Range  for  Electrons  in  Aluminum.  Low  Energy  Values  (Blue  Diamonds)  from 
Ashley  et  al.  (Reference  17).  High  Energy  Values  (Red  Squares)  from  NIST  ESTAR*  Database. 
Dashed  Line  is  Bi-exponential  Fit. 

For  materials  where  no  suitable  data  is  available,  a  mono-exponential  form  can  be  generated 
using  Feldman’s  empirical  relationships16,  connecting  b  and  n  to  atomic  data. 

b  =  250w/pmassZcl/2  (23) 

q  =  1  -2/(1  -  0.29  log10  Z)  (24) 

where  W  is  the  atomic  or  molecular  weight  of  the  material,  Z  is  the  atomic  number,  and  pmass  is 
the  mass  density  in  gm  cm'  .  Then  E  is  in  keV  and  R  is  in  angstroms.  The  stopping  power  is  then 


National  Institute  of  Standards  and  Technology  Electron  STopping  power  And  Range  database  available  at 
http://www.nist.gov/pml/data/star/index.cfm 


12 

Approved  for  public  release;  distribution  is  unlimited. 


obtained  indirectly  with  the  above  equation.  Theoretical  estimates  of  the  stopping  power  for  a 
number  of  materials  are  available  from  Ashley  et  al.  Comparison  of  these  values  with  those 
implied  by  the  range  data  showed  significant  discrepancies,  particularly  for  those  materials  fit 
using  Feldman’s  fonnula.  The  best  approach  is  to  fit  the  four  parameters  in  the  equation  for  R 
directly  to  the  stopping  power  data. 

S  =  (q1b1Eqi_1 +q2b2Eq2_1)  1  (25) 


The  secondary  electron  yield  depends  on  angular  distribution  of  the  incident  electrons.  For  low 
energy  electrons,  all  the  energy  is  deposited  near  the  surface  regardless  of  incident  angle,  so  the 
yield  is  independent  of  incident  angle.  For  high  energy  electrons  the  stopping  power  is  constant 
in  the  near-surface  region,  so  the  yield  is  proportional  to  the  path  length  near  the  surface,  or 
inversely  proportional  to  cos \|/.  It  follows  that  the  yield  for  high  energy  isotropic  electrons  is 
double  that  for  normally  incident  electrons  of  the  same  energy.  The  isotropic  yield  curve  makes  a 
transition  from  being  equal  to  the  normal  incidence  curve  at  low  energy  to  double  at  high  energy. 
The  maximum  yield  for  isotropic  electrons  is  somewhat  higher  than  for  normal  electrons,  and 
occurs  at  somewhat  higher  energy. 

2.2.2  Secondary  Electron  Emission  Due  to  Ion  Impact 

Secondary  emission  of  electrons  due  to  ion  impact  can  be  treated  in  a  way  similar  to  that  for 
electron  impact.  The  yield  is  given  by 


Y: 


le 


(26) 


The  stopping  power  is  assumed  to  be  independent  of  path  length  x  over  the  thickness,  d,  of  the 
sample.  Above  10  keV,  the  stopping  power  is  approximately  proportional  to  the  velocity,  and  at 
higher  energies  (above  -200  keV)  it  is  inversely  proportional  to  the  velocity  .  Below  10  keV,  it 
is  lower. 


dE 


dx 


P(e)e 


1/2 


1  +  E/E 


(27) 


max 


Emax  is  the  energy  at  the  maximum  in  the  yield  curve.  This  is  approximately  100  keV  for  most 
materials. 


f  0.25 


P(E)  = 


E  <  0.476  keV 

0.476  keV  <  E  <  10  keV 
10  keV  <  E 


(28) 


As  the  CSDA  range  is  at  least  comparable  to  the  mean  depth  of  secondary  electron  production 
for  protons  above  1  keV  (140  angstroms  of  aluminum),  we  consider  the  proton  yield  to  vary 
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inversely  with  the  cosine  of  the  incident  angle,  so  that  the  isotropic  incidence  yield  is  double  the 
normal  incidence  yield.  This  is  usually  an  adequate  treatment,  because  the  ion  fluxes  are  almost 
always  low,  and  the  effect  of  secondary  emission  is  to  mildly  augment  the  incident  ion  flux, 
rather  than  canceling  it  as  is  the  case  with  incident  electrons. 

The  secondary  emission  properties  due  to  the  impact  of  ions  other  than  protons  are  assumed  to  be 
identical  to  the  proton  values  for  the  same  energy. 

The  secondary  electron  yield  curve  for  protons  incident  on  aluminum  is  shown  in  Figure  6.  As 
can  be  seen  by  comparison  with  Figure  7,  the  yield  closely  follows  the  stopping  power. 


- 1.36  E1/2/(1+E/40) 

□  Cousinie  et  al. 

[1959] 

*  Hill  etal.  [1939] 

_ 


Figure  6.  Secondary  Electron  Emission  by  Aluminum  for  Proton  Impact  at  Normal  Incidence; 
Experimental  Points  as  Indicated  in  the  References. 19, 20’ 21, 22 
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Figure  7.  Stopping  Power  for  Protons  on  Aluminum.  (From  NIST  PSTAR  Database.1) 

2.2.3  Backscatter 

Backscattered  electrons  are  those  emitted  from  the  surface  with  energies  above  50  eV.  Their 
energy  distribution  is  usually  peaked  close  to  the  primary  incident  energy  and  they  may  be 
considered  as  reflected  electrons. 

23  24  25 

Nascap-2k  uses  a  backscattering  theory  based  on  that  of  Everhart"  as  extended  by  McAfee  . 
It  assumes  a  single  scattering  in  accordance  with  the  Rutherford  cross-section  and  the  Thomson- 
Widdington  slowing  down  law, 


dE 

dx 


oc  E 


(29) 


(valid  for  most  metals  for  E  >  10  keV).  For  normal  incidence  the  backscattering  coefficient  is 
given  by 


0  =  1 


(30) 


where  a  is  taken  to  be  0.0375  Z  and  where  Z  is  the  atomic  number  of  the  material.  This 
expression  matches  the  experimental  data. 

The  large-angle  scattering  theory,  together  with  Monte  Carlo  data  and  experiments  by  Darlington 
and  Cosslett,  indicate  that  the  angular  dependence  of  backscattering  is  well-described  by 

0(¥)  =  0oexP(rli(l-cos\)/))  (31) 


1  National  Institute  of  Standards  and  Technology  Proton  STopping  power  And  Range  database  available  at 
http://www.nist.gov/pml/data/star/index.cfm 
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where  the  value  of  r)j  is,  within  the  uncertainty  in  the  data,  what  would  be  obtained  by  assuming 
total  backscattering  at  glancing  incidence,  r^  =  -log  r|o.  The  net  albedo  for  an  isotropic  flux  is 
then 


A  _21~ho  (1-log  ho) 

(!°gho)2  (32) 

97 

As  the  energy  is  decreased  below  10  keV,  the  backscattering  increases.  Data  cited  by  Shimizu" 
indicate  an  increase  of  about  0.1,  almost  independent  of  Z.  This  component  of  backscattering  can 
be  approximated  by 


5r|0  =0.1exp[-E/5  keV] 

At  very  low  energies,  the  backscattering  coefficient  becomes  very  small  and,  below  50  eV, 
becomes  zero  by  definition,  as  all  emitted  electrons  are  considered  secondaries.  This  can  be 
taken  into  account  by  a  factor  of 


H(E  -  50eV) 


log  20 


log 


50eV 


(34) 


where  H  is  the  Heaviside  step  function.  The  formula  for  energy-dependent  backscattering, 
incorporating  these  assumptions,  is  then 


ho 


H(1-E)H(E  -  0.05) 


f  1  A 

log  20 


log 


0.05 


+  H(E-1) 


3  1  e-E/5 


10 


+  1- 


\0.037Z  ) 


v  ey 


(35) 


where  energies  are  measured  in  keV. 

2.2.4  Photoemission 


Photoelectrons  are  those  ejected  from  the  surface  due  to  the  solar  ultraviolet  radiation.  Usually 
the  quantity  known  is  the  yield,  or  number  of  electrons  emitted  for  a  surface  normally  exposed  to 
the  solar  spectrum,  at  an  “earth  distance”  from  the  sun  (Ysun).  The  photocurrent  from  a  surface 
exposed  to  the  sun  at  an  angle  \|/sun  is  given  by  the  fonnula 

i photo  =  (Area  exposcd)Ysun  cos ysun  (36) 


This  assumes  that  the  yield  per  photon  is,  on  average,  independent  of  \|/sun. 

A  surface  element  is  taken  to  be  sunlit  if  its  centroid  is  sunlit,  i.e.,  if  a  vector  from  the  centroid  in 
the  sun  direction  does  not  intersect  any  other  surface  element. 
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2.3  Surface  Charging  from  Tracked  Currents  and/or  Current  Balance 

In  short  Debye  length  environments,  such  as  encountered  in  low-Earth  orbit,  the  orbit-limited 
approximation  greatly  overestimates  currents  of  attracted  particles.  This  occurs  because 
potentials  that  decay  faster  than  inverse  square  (sheath  potentials)  lead  to  an  effective  potential 
barrier  to  attracted  charged  particles  when  angular  momentum  is  taken  into  account.  Nascap-2k 
can  compute  the  charging  of  spacecraft  surfaces  using  tracked  currents,  either  by  themselves  or 
in  conjunction  with  analytically  computed  currents.  Calculations  of  typical  interest  include: 

1.  Floating  potential  and/or  differential  charging  of  a  negative  spacecraft,  whose  ion  current 
may  be  enhanced  and  asymmetrized  due  to  ram/wake  effects. 

2.  Floating  potential  of  a  biased  spacecraft  that  has  both  ion  and  electron  sheaths.  The 
electron  current  may  be  additionally  limited  by  magnetic  field  effects. 

3.  Dynamic  potential  of  a  spacecraft  with  time-dependent  bias. 

In  “Surface  Charging”  calculations  with  “Tracked  Ion  and  Analytic  Electron  Currents” 
(“Auroral”  environment  only),  the  analytic  electrons  used  are  high  energy  electrons  computed 
from  the  Fontheim  formula.  In  “Time  Dependent  Plasma”  calculations,  the  analytic  electrons 
used  are  computed  by  the  convected  Maxwellian  formula. 

2.4  Potentials  in  Space 

Nascap-2k  uses  a  finite  element  method  to  solve  Poisson’s  equation  -  s0V2(|)  =  p(x)  for 
electrostatic  potentials  throughout  space.  Nascap-2k  uses  the  variational  fonn 


o  =  A{ 

• 

dV 

6<t)  IJ 

L2  8oJ 

+  {<j)V<|)*dS 


(37) 


The  first  tenn  in  the  volume  integrand  corresponds  to  the  Laplacian  operator;  the  second  tenn  is 
the  space  charge  contribution.  As  surface  potentials  are  held  fixed  when  solving  for  potentials  in 
space,  the  surface  tenn  appearing  in  Equation  (37)  is  not  used  in  the  potential  solution,  but  may 
be  used  later  to  determine  the  normal  component  of  electric  field  (and  thus  the  surface  charge) 
for  each  surface  element. 

The  solution  is  subject  to  fixed  potential  boundary  condition  on  the  spacecraft  surfaces  and  either 
zero  potential  or  monopole  potentials  on  the  grid  boundary. 

Available  Space  Charge-Density  Models 

A  number  of  models  of  the  charge  density  are  available  to  study  different  phenomena. 

In  the  following,  the  symbol  A,D  is  the  plasma  Debye  length  and  g  is  the  plasma  density  reduction 
factor  computed  by  geometric  wake  model  (0<g<l). 

Laplace.  The  Laplacian  space  charge  option  specifies  that  the  charge  density  is  zero  and  is 
appropriate  for  very  low  density  plasmas. 
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(38) 


i.e.,  charge  exists  only  on  object  surfaces  and  external  boundaries,  as  detennined  by  the 
boundary  conditions.  “Space  charge”  iterations  may  still  be  required,  however,  due  to  the 
treatment  of  surface  electric  fields. 


Linear  (Debye  Shielding).  The  Linear  space  charge  option  solves  the  Helmholtz  or  Debye  - 
Huckel  equation.  This  model  is  rarely  applicable  to  spacecraft  because  it  assumes  all  potentials 
are  much  less  than  the  plasma  temperature,  which  is  only  true  in  low  density  plasma  where  the 
laplace  approximation  is  valid. 


(39) 


Nonlinear.  The  Nonlinear  space  charge  model  is  appropriate  for  most  low-Earth-orbit  type 
plasmas.  It  accounts  for  space  charge  acceleration  and  convergence  in  a  manner  based  on 
spherical  collection  (Langmuir-Blodgett  problem).  Poisson’s  equation  is  solved  with  space 
charge  given  by: 


-V2(j)  = 


4) 

(  max(l,  c((j), 

E 

^o/g 

vi+V4^/e3/2  J 

(40) 


This  model  smoothly  interpolates  between  Debye  screening  at  low  potential  and  an  accelerated 
distribution  with  particle  convergence  at  high  potentials.  The  convergence  factor  is  computed  in 
terms  of  local  information  and  problem  parameters.  The  convergence  model  was  developed  by 
numerically  solving  the  Langmuir-Blodgett  problem  for  collection  by  a  high-voltage  sphere  and 
fitting  the  result  to  an  analytic  fonn.  An  excellent  fit  was  found  with 

{RSh  A)2  =  2.29(|e|  A,d/0^1'262 \Oj  0p°'509 

C(<(),  |E|)  =  min((Rsh  /r  f  ,3 -545|c|>/e| 3/2 ) 
where  E  is  the  local  electric  field. 


(41) 

(42) 


Frozen  Ion.  The  Frozen  Ion  formulation  is  intended  for  short  timescale  (typically  sub¬ 
microsecond)  problems  for  which  it  is  a  good  approximation  to  assume  that  ions  remain 
stationary  and  at  ambient  density  (“ion  matrix”  approximation),  but  electrons  achieve  barometric 
equilibrium.  This  approximation  is  intended  primarily  for  negative  potentials,  for  which  the 
charge  density  is  zero  at  zero  potential  and  for  which  at  negative  potential  (electrons  repelled)  the 
charge  density  becomes  positive  and  asymptotically  reaches  the  ion  density.  For  computational 
convenience  we  extend  the  charge  density  linearly  to  positive  potentials,  as  an  exponential 
increase  in  electron  charge  would  be  unphysical. 
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(43) 


P  =  gne 
P  =  - 


())  <  0 
(j)  >  0 


Considerations  for  implementing  this  formulation  stably  on  a  coarse  mesh  are  discussed  in 
Section  3.6.5. 


Barometric.  The  barometric  space  charge  model  is  for  cases  in  which  all  the  surfaces  are  at 
potentials  comparable  to  or  below  the  plasma  temperature  and  there  is  a  region  of  low  density, 
such  as  a  plasma  wake.  The  total  charge  density  is  the  sum  of  the  ion  and  electron  charge 
densities. 


p  =  pe  +Pi 


(44) 


The  ion  density  is  the  plasma  density  decreased  by  the  wake  factor. 

Pi  =  g  en 

The  electron  density  is  expanded  about  barometric  equilibrium: 


(45) 


Pe  =-pi< 


1  +  (<t>  -  <t>b  )/0  4>>4> 

exp(4-(j)b)/0)  <|><<|> 
<l>b  =  01n(g) 


(46) 

(47) 


In  a  dense,  short  Debye  length  plasma,  the  requirement  that  the  ion  and  electron  densities  be 
nearly  equal  gives  strictly  barometric  potentials.  In  plasmas  with  a  longer  Debye  length  strict 
quasineutrality  does  not  hold,  so  that,  for  example,  the  potential  in  a  wake  is  considerably  less 
negative  than  barometric,  and  the  wake  is  therefore  electron  rich. 

Note  that  this  model  differs  from  the  “Nonlinear”  model  in  that  it  lacks  an  approximation  to 
attracted  species  density  in  sheaths  and  does  not  allow  for  ion  depletion  in  positive  potential 
regions. 

Full  Trajectory  Ions.  Ion  densities  are  calculated  from  steady-state  ion  trajectories.  Electrons 
are  barometric.  This  approximation  is  used  to  model  charge  density  for  a  moving  spacecraft 
when  current  collection  in  the  wake  is  important  and  the  neutral  wake  approximation  is 
inadequate — such  as  when  high  potential  surfaces  are  in  the  wake. 


—  =  —  (i-exP((4>-4>b)/e)) 


s0  e 

4>b  =  ein| 


leny 


(48) 


Here,  p;  is  calculated  by  tracking  ions. 
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Plume  Ion  Density.  Ion  densities  are  initially  computed  from  an  imported  plume  map  file. 
Electrons  are  barometric.  The  quantity  p/s0  is  computed  using  the  same  fonnula  as  for  full 

trajectory  ions,  where  pi  is  computed  by  summing  the  contributions  from  the  thruster  plumes 
instead  of  computed  by  tracking  macroparticles.  After  the  initial  iteration,  pi  is  computed  by 
summing  the  main  beam  contribution  from  the  thruster  plumes  and  the  charge  exchange 
contribution  from  generating  and  tracking  charge  exchange  ions. 

Hybrid  PIC.  This  algorithm  is  used  for  timescales  (typically  sub-millisecond)  on  which  it  is 
practical  to  treat  ion  motion,  but  electrons  may  be  considered  in  barometric  equilibrium.  The 
total  charge  density  is  the  sum  of  the  ion  and  electron  charge  densities. 

P  =  Pe+Pi  (49) 

The  ion  density  is  computed  from  ion  macroparticles  that  move  in  the  local  electric  fields  at  each 
timestep.  The  electron  charge  density  is  given  by  the  same  fonnula  as  for  the  Barometric  model 
above  with  the  barometric  potential  given  by 


c|)b  =  01n 


fen) 


(50) 


Full  PIC.  For  this  option,  it  is  assumed  that  both  the  electron  and  ion-charge  densities  were 
stored  during  particle  tracking. 


_  _  tracked  /  c  i  \ 

P  =  P  (51) 

2.5  Macroparticle  Creation  and  Tracking 

Particle  tracking  is  used  to  study  sheath  currents,  to  study  detector  response,  to  generate  steady- 
state  charge  densities,  and  to  generate  space  charge  evolution  for  dynamic  calculations.  The 
investigation  of  detector  response  is  applicable  in  both  tenuous  and  dense  plasmas.  However,  the 
other  phenomena  are  only  applicable  in  plasmas  with  Debye  lengths  no  more  than  the  spacecraft 
size. 

In  a  dense  plasma  with  high  potentials,  the  disturbed  region  is  referred  to  as  a  sheath.  The 
“sheath  surface”  represents  a  sharp  demarcation  between  a  low  potential  exterior  region 
containing  neutral  “undisturbed  plasma,”  and  a  high  potential  interior  region  from  which  one 
species  (ions  for  positive  potentials;  electrons  for  negative  potentials)  is  excluded. 

In  a  dense  plasma,  the  total  current  and  its  distribution  over  the  spacecraft  surface  depend  on  the 
plasma  environment,  the  surface  potentials,  the  sheath  structure,  the  spacecraft  velocity,  and 
ambient  and  spacecraft  generated  magnetic  fields.  For  a  large  object  (many  Debye  lengths  in 
size)  at  low  potential  (at  most,  a  few  times  the  plasma  temperature),  Debye  screening  minimizes 
the  extent  of  the  object’s  influence  on  the  plasma.  Because  the  sheath  is  “thin,”  the  sheath  surface 
lies  near  the  spacecraft  surface.  The  current  attracted  to  each  surface  element  is  then  given  by  its 
area  times  the  plasma  thennal  current  of  the  appropriate  species.  For  an  object  (not  necessarily 
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large)  at  high  potential,  however,  the  extent  of  perturbed  plasma  is  far  greater.  A  first  estimate  of 
the  sheath  thickness  is  given  by  the  Child-Langmuir29,30  distance. 


D 


CL 


8s0  frc 

9en  V  0 


(52) 


The  plasma  thennal  current  then  flows  through  an  area  much  larger  than  the  spacecraft  surface. 
The  current  may  be  focused  by  the  fields  within  the  sheath,  which  affects  the  distribution  of 
current  over  the  spacecraft  surface.  This  focusing  causes  the  increase  in  space  charge  density  due 
to  particle  acceleration  modeled  in  the  “Nonlinear”  space  charge  density  model. 

Macroparticles  are  generated  in  the  Particle  Generation  step  of  the  script  and  tracked  in  the 
Track  step  of  the  script.  Macroparticles  can  be  generated  at  a  well-defined  sheath  surface,  at  the 
boundary  of  the  computational  space,  throughout  the  volume  of  space,  at  surface  elements,  or  at 
user  specified  locations.  As  macroparticles  are  tracked,  their  charge  density  in  the  volume 
elements  they  pass  through  may  or  may  not  be  retained,  depending  on  user  choices.  The  current 
to  any  surface  element  hit  by  the  macroparticles  is  retained. 

2.5.1  Generation  of  Macroparticles  at  Sheath  Surface 

The  attracted  species  diffuses  across  the  sheath  surface  at  a  rate  given  by  the  plasma  thermal 
current. 


/  e0 

J‘=en&  (53 

The  current  density  to  the  sheath  edge  is  computed  by  assuming  the  sheath  to  be  a  perfectly 
absorbing  spherical  surface. 

The  charged  particles  constituting  this  current  undergo  motion  consistent  with  the  electric  and 
magnetic  fields  inside  the  sheath.  The  sheath  surface  is  taken  as  the  equipotential  surface  at 
(j)  =  ±01n2.  This  choice  is  made  because  the  attracted  species  is  absorbed  by  the  sheath,  so  we 
have  only  the  inward  moving  component,  comprising  half  the  ambient  density.  The  repelled 
species,  whose  density  satisfies  n  =  na  exp(-  |(J>/0|),  must  also  be  at  half  the  ambient  density, 
leading  immediately  to  the  sheath  potential. 

Macroparticles  representing  this  sheath  current  can  be  generated  and  tracked  to  detennine  the 
current  to  the  spacecraft  and  its  distribution.  The  initial  velocity  is  the  average  velocity  of  a 
charged  particle  crossing  the  sheath  boundary.  In  the  absence  of  spacecraft  motion  or  a  magnetic 

|2e0 

field,  this  is  v  =  J -  in  the  direction  of  the  local  electric  field.  More  precisely,  for  a  sheath 

V  Ttm 

segment  on  the  ram  side 
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where  the  plus  sign  is  for  ions  and  the  minus  for  electrons  and  vx  is  the  component  of  the 
spacecraft  velocity  along  the  local  electric  field.  In  the  absence  of  spacecraft  motion  and  for  a 
sheath  segment  in  the  wake, 


With  a  unifonn  magnetic  field 


where  the  drift  velocity  tenn  vdrift 
square  root  negative. 


E  x  B 

— —  is  not  included  if  it  would  make  the  argument  of  the 


This  model  is  appropriate  for  high  potentials.  When  either  the  thennal  distribution  of  particle 
velocities  or  the  spacecraft  velocity  is  important,  macroparticles  should  be  tracked  from  the 
boundary  of  the  computational  space. 


2.5.2  Generation  of  a  Thermal  Distribution  at  the  Computational  Space  Boundary 


There  are  two  models  for  creation  of  macroparticles  at  the  boundary  of  the  computational  grid.  In 
both  models,  the  plasma  is  assumed  to  be  Maxwellian.  The  two  models  are  inconsistent  and 
should  not  be  used  together. 


Boundary  Injection 


In  the  simplest  model,  macroparticles  are  created  with  charge  equal  to  the  plasma  thennal  current 


times  the  area  times  the  time  interval,  q  =  jthAAt,  and  velocity  equal  to  J -  normal  to  the 

V  7im 


boundary,  so  that  they  represent  a  density  of  n/2.  These  can  each  be  split  into  eight  outwardly 
moving  macroparticles,  to  approximate  a  thennal  distribution.  The  split  macroparticles  have  an 


additional  velocity  of  ±0.707 


along  each  axis  of  a  randomly  generated  coordinate  system. 


The  splitting  is  done  in  the  plasma  frame  of  reference  in  order  to  simulate  the  correct  momentum 
and  energy  distribution  for  a  drifting  Maxwellian  when  transfonned  back  to  the  spacecraft 
reference  frame. 
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Boundary 

In  this  model,  macroparticles  that  represent  a  specified  fraction  of  the  thennal  distribution  in 
each  of  the  three  spatial  directions  are  created.  At  each  boundary  emission  point,  several 
macroparticles  are  created  that  sample  the  velocity  distribution  function  in  each  of  the  three 
spatial  directions.  The  user  specifies  which  fraction  of  the  distribution  function  each 
macroparticle  is  to  represent  and  the  code  computes  the  range  of  velocities  (vf’2y,z )  which  would 
give  the  specified  fractions,  ip. 


v2  v2  v2  v2  V2  V2 
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Several  macroparticles  are  created,  each  with  velocity  the  sum  of  the  average  velocity  of  its 
fraction  of  the  distribution  and  the  ram  velocity  (negative  of  the  spacecraft  velocity). 
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2.5.3  Generation  of  Macroparticles  Throughout  the  Computational  Space 
Initialization 


Macroparticles  can  be  created  throughout  the  computational  space  to  represent  a  unifonn 
distribution  of  charge.  The  weight  of  these  macroparticles  is  decreased  by  any  geometric  wake 
factor.  Each  macroparticle  can  be  split  into  eight  outwardly  moving  macroparticles  so  as  to 
approximate  a  thennal  distribution.  The  velocities  have  components  of  ±0.707  ^eO/m  along 
each  axis  in  a  randomly  oriented  coordinate  system.  The  splitting  is  done  in  the  plasma  frame  of 
reference  in  order  to  simulate  the  correct  momentum  and  energy  distribution  of  a  drifting 
Maxwellian  when  transfonned  back  to  the  spacecraft  reference  frame. 

Charge  Exchange 

Macroparticles  can  also  be  created  throughout  the  computational  space  to  represent  charge 
exchange  current.  Charge  exchange  ions  are  those  created  by  a  charge  transfer  collision  between 
a  fast  ion  and  a  slow  neutral  atom,  such  as  occurs  in  ion  thruster  plumes.  The  current  per  unit 
volume  associated  with  each  macroparticle  is  the  product  of  the  charge  exchange  cross  section, 
the  neutral  density,  and  the  main  beam  ion  current  density. 

Jcex  C  C5cex  ^neutral  Abeam  Vbeam-  (59 


Each  component  of  the  initial  velocity  is  given  by 


2e0.3448 


m 


where  0.3448  is  4000  K  in 


electron  volts  and  C,  is  a  random  number  between  -1  and  1. 
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The  neutral  density  is  the  sum  of  the  un-ionized  propellant  from  the  thrusters,  the  gas  flowing 
through  the  neutralizers,  and  the  background  gas. 


n  unionized  gas 

neutral  thrusters  neutralizers 


neutral 


(60) 


The  contribution  from  each  thruster  is  given  by  the  solid  angle  subtended  by  the  thruster  grid,  the 
flow  rate  and  the  temperature. 
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The  contribution  from  each  neutralizer  is  given  by  the  flow  rate,  the  temperature,  the  angle 
between  the  neutralizer  axis  and  the  line  of  sight  from  the  neutralizer,  and  the  distance  from  the 
neutralizer. 
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2.5.4  Generation  of  Macroparticles  at  Surface  Elements 

Macroparticles  can  be  created  on  surface  elements  to  determine  if  emitted  current  will  return  to  a 
surface  of  the  object  or  leave  the  computational  space,  or  to  determine  the  current  incident  to  a 
detector  from  a  convected  Maxwellian  plasma. 

Emitted  Current 

Macroparticles  that  represent  a  user  specified  emitted  current  and  range  of  angles  and  kinetic 
energies  are  created  on  surface  elements.  The  macroparticles  are  emitted  with  a  flat  distribution 
in  velocity,  azimuthal  angle,  and  cosine  of  the  polar  angle.  A  full  distribution  of  macroparticles  is 
created  at  a  user  specified  number  of  points  on  each  emitter  surface  element.  The  emitted  current 
is  subtracted  from  the  total  current  to  the  emitting  surface  element. 

Current  to  Detectors  (Reverse  Trajectory  Tracking) 

Another  reason  to  create  macroparticles  on  surface  elements  is  to  compute  the  current  to  those 
surface  elements  from  a  flowing  thennal  plasma.  This  technique  is  applicable  when  the  current 
of  interest  is  composed  of  charged  particles  with  velocity  (direction,  magnitude,  or  both)  far  from 
the  mean  of  the  distribution  function,  the  sheath  is  very  large,  or  there  is  some  other  condition 
such  that  tracking  particles  inward  from  a  sheath  boundary  does  not  work  well.  Under  these 
conditions,  reverse  tracking  can  be  used  to  detennine  the  original  velocity  of  the  incident 
particles,  or  if  they  are  blocked.  The  current  is  then  given  by  an  integral  over  the  environment 
distribution  function  of  the  incident  species.31 

Macroparticles  that  represent  a  user  specified  field  of  view  are  created  at  a  user  specified  number 
of  points  on  each  detector  surface  element.  The  user-specified  number  of  macroparticles  are 
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distributed  evenly  in  velocity,  azimuthal  angle,  and  cosine  of  the  polar  angle.  Each  macroparticle 
has  a  weight  appropriate  to  its  fraction  of  the  total  incident  current. 

When  these  macroparticles  exit  the  computational  space,  the  plasma  properties  and  each 
macroparticle’s  weight  and  final  velocity  are  used  to  assign  an  incident  current  to  the  originating 
surface  element. 

2.5.5  Particle  Tracking 

Macroparticles  are  tracked  until  they  leave  the  computational  grid,  hit  a  surface,  exceed  the 
number  of  allowed  steps,  or  exceed  the  computational  tracking  time.  They  are  tracked  in  the 
local  electric  field  and  the  magnetic  field  (ambient  unifonn  and  spacecraft-generated  dipoles) 
using  a  third  order  energy  conserving  algorithm. 

As  the  macroparticles  pass  through  volume  elements,  their  charge  and/or  current  can  be 
optionally  saved  to  volume  elements  and  nodes.  This  charge  and/or  current  can  subsequently  be 
used  in  the  space  charge  fonnulas  to  compute  potentials  in  space. 

2.6  Geometric  Wake 

The  exclusion  of  plasma  in  the  wake  of  a  rapidly  moving  spacecraft  is  also  important.  That 
portion  of  the  sheath  surface  that  lies  in  the  wake  region  has  a  reduced  current  passing  through  it. 
The  volume  of  space  in  the  wake  has  a  reduced  charge  density.  Conversely,  on  the  ram  side  of 
the  spacecraft,  there  is  enhanced  current  and  charge  density. 

A  spacecraft  creates  a  wake  when  its  velocity  is  comparable  to,  or  larger  than,  the  thermal 
velocity  of  the  ambient  ions.  In  low  Earth  orbit,  the  spacecraft  velocity  is  7500  ms'1.  The  thennal 
velocity  of  a  0.3  eV  oxygen  ion  is  only  about  2000  ms'1.  Thus  the  spacecraft  travels  a  few  of  its 
own  radii  before  ions  can  fill  in  behind  it.  The  electrons  fill  in  more  rapidly.  However,  the 
density  to  which  electrons  can  accumulate  is  limited  by  the  space  charge  of  the  electrons  already 
in  the  wake.  Thus,  except  in  regions  where  the  density  is  extremely  low,  the  electron  density  is 
only  very  slightly  greater  than  the  ion  density. 

For  each  point  in  space,  the  wake  module  calculates  what  fraction  of  the  thennal  distribution  of 
velocities,  (in  the  spacecraft  frame,  a  Maxwellian  displaced  by  the  spacecraft  velocity  vector)  is 
not  blocked  by  the  spacecraft. 

The  neutral  wake  density  is  given  by 


M=  kUfZ)!  f0(v,nydv  \dD. 


(63) 


where  f0(v,f2)  is  the  unperturbed  distribution  function  for  a  drifting  Maxwellian,  and  g(x,Q)  has 
value  “0”  if  a  ray  starting  from  x  and  going  in  the  direction  Q  would  strike  the  spacecraft  and  “1” 
if  it  would  not. 
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2.7  Plumes 


The  ion  density  and  velocity  due  to  an  ion  thruster  main  beam  can  be  read  into  Nascap-2k  and 
displayed.  Multiple  thrusters  (each  with  its  own  location  and  axis  direction)  can  be  specified, 
although  at  present  all  thrusters  have  the  same  type  of  plume.  Neutralizers  can  be  specified  as 
well.  The  plume  currents  are  used  in  the  computation  of  the  potentials  in  space  (“Plume  Ion 
Density”  space  charge  model)  and  charge  exchange  currents. 

2.8  Transverse  Surface  Currents 


The  transverse  surface  current  algorithm  in  Nascap-2k  is  used  to  calculate  the  current  flowing 
along  an  antenna  element  or  other  object  surface  in  response  to  time-varying  applied  potentials 
and/or  active  or  passive  current  sources.  The  model  accounts  correctly  for  both  local  capacitance 
and  incident  plasma  current.  The  change  in  charge  on  a  surface  element  during  a  timestep  is  that 
which  is  needed  to  accomplish  the  change  in  surface  electric  field,  which  is  obtained  from  the 
Nascap-2k  potential  calculation.  The  current  to  each  surface  element  consists  of  the  transverse 
surface  current  from  neighboring  surface  elements  and  the  current  provided  by  the  plasma,  which 
is  provided  by  the  Nascap-2k  PIC  simulation.  The  current  continuity  equation  is  solved  using  a 
pseudopotential  approach.  As  a  boundary  condition,  one  surface  element  of  each  conductive 
element  (i.e.,  surfaces  associated  with  a  given  Nascap-2k  conductor)  is  specified  as  connected  to 
the  biasing  power  supply,  and  the  solution  provides  the  required  current. 


The  basic  equation  to  be  solved  is  V  •  J  +  -A-  =  0 ,  where  J  is  the  surface  current  density,  and  p  is 

the  surface  charge  density.  Note  that  the  solution  for  the  current  is  non-unique  to  the  addition  of 
any  divergence-free  current  field.  We  assume  that,  provided  appropriate  boundary  conditions  are 
implemented,  such  circulating  currents  are  not  of  concern  to  the  problem  being  solved. 


2.9  Propagating  Fields 

The  transverse  surface  currents,  the  volume  ion  currents  computed  during  particle  tracking,  and 
the  volume  electron  currents  computed  from  tracking  or  saved  in  the  database  by  an  external 
code,  are  a  source  of  propagating  electromagnetic  fields.  Nascap-2k  can  compute  the  magnetic 
field,  the  vector  potential,  the  transverse  electric  field  (rate  of  change  of  the  vector  potential),  and 
the  Poynting  vector  from  these  currents. 


.  trans 
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p(r)  =  £oC2^xB(r) 
dt 


(66) 
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where  the  surface  integrals  are  over  the  spacecraft  surfaces  and  the  volume  integrals  are  over  all 
space  outside  the  spacecraft. 


2.10  Electron  Currents  (Using  Pseudopotential) 

Nascap-2k  includes  an  external  code  that  provides  a  non-PIC  method  to  calculate  volume 
electron  currents  that  satisfy  the  equation  of  continuity,  along  with  other  physical  requirements, 
and  are  consistent  with  electron  densities  computed  by  an  ion  dynamics  Nascap-2k  simulation. 
The  method  can  be  used  to  calculate  electron  currents  within  and  near  the  sheath  about  a  VLF 
(Very  Low  Frequency)  antenna  or  other  high-voltage  object  at  frequencies  that  are  low  compared 
with  both  the  electron  plasma  frequency,  wpe,  and  the  electron  cyclotron  frequency,  cocc.  The 
currents  can  then  be  used  to  calculate  electromagnetic  radiation,  ohmic  heating,  etc. 


This  external  code  uses  a  pseudopotential  approach  to  the  computation  of  volume  electron 
currents.  Local  equilibrium  electron  densities  are  generated  for  each  volume  element  as  part  of 
the  Hybrid  PIC  ion  dynamics  simulation,  or,  if  not  available,  may  be  calculated  by  the  module  as 


1  +  ^/  cp  >  0 
e1^1'  (p<0 


(67) 


Time  derivatives  of  the  electron  density  are  the  main  drivers  of  volume  electron  currents.  Space 
outside  the  calculation  boundary  can  act  as  either  a  source  or  a  sink  for  electrons,  and  object 
surfaces  may  act  as  electron  sinks. 


The  electron  current,  jeiec,  must  satisfy 


V  ’  jeiec  + 


dpelec 

dt 


=  0  .  To  solve  this  we  assume  that  jeiec  is 


proportional  to  the  gradient  of  a  “pseudopotential,”  vp:  jelec  =  oVv|/ ,  where  the  conductivity 

tensor,  o,  depends  on  the  electron  density  and  magnetic  field.  Note  that  the  solution  for  the 
current  is  non-unique  to  the  addition  of  any  divergence-free  current  field.  The  above  equations 
can  be  combined  as 


dp 


elec 


-V  -aV\|/  = 

dt  .  (68) 

By  standard  finite  element  treatment,  this  equation  is  equivalent  to  minimizing  the  function 


dV 


— Vv|/  -(tV\|/  -\\ ) 


dPelec 

dt 


(69) 


The  conductivity  tensor  to  be  applied  to  the  pseudopotentials,  o,  is  assumed  to  be  the  same  as  the 
nonnal  conductivity  tensor,  which  for  B||z  is  given  by 


27 

Approved  for  public  release;  distribution  is  unlimited. 


°std  - 


1  +  (coct)2 


1 

*Y 

0 


COcT 


-COcT  0  N 

1  0 

0  l  +  (cocx)2y 


(70) 
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We  take  a0  = - and  —  =  — — ,  so  that  cocx  =  — - .  The  otf-diagonal  (Hall)  terms  are  omitted 

m  x  4  cone 

pc 

for  the  pseudopotential  solution  and  restored  when  calculating  the  currents.  To  obtain  the 
conductivity  tensor  for  an  arbitrary  magnetic  field  we  apply  the  transfonnation 


°  =  TTostdT 

where  the  transformation  T  is  given  by 
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with  cos  0  =  -j — r ,  tan  m  =  — . 
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3  NUMERIC  MODELS 


This  section  describes  the  numeric  techniques  Nascap-2k  uses  to  implement  the  physics  models 

described  in  the  previous  section. 

•  The  orbit  limited  currents  are  computed  by  integrating  over  the  distribution  function.  An 
implicitized  version  of  the  Boundary  Element  Method  is  used  to  compute  the  evolution  of 
surface  potentials. 

•  The  finite  element  approach  is  used  to  solve  Poisson’s  equation  for  the  spatial  potential. 
Implicitization  and  charge  stabilization  are  used  for  stability.  The  finite  element  grid  uses 
strictly  continuous  electric  field  interpolants  for  accurate  particle  tracking. 

•  Macroparticle  tracking  is  used  for  computing  surface  and  volume  currents  and  volume 
densities.  Macroparticle  splitting  is  used  for  improved  accuracy. 

3.1  Environment  Integrals 

In  section  2.1  we  noted  that  the  analytic  current  density  to  a  surface  at  potential  (j)  is  given  by 
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where  L  is  0  for  the  repelled  species  and  |(j)|  for  the  attracted  species,  the  “+”  sign  is  for  ions  and 
the  sign  is  for  electrons,  f(E)  is  the  distribution  function  at  infinity,  and  F(E)  is  the  flux  at 
infinity. 

Below  we  describe  how  these  integrals  are  done. 

The  environment  parameters  used  are  those  valid  at  the  beginning  of  the  timestep. 

3.1.1  Maxwellian  Currents 
The  Maxwellian  integrals  are  simple. 
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which  for  the  repelled  species  simplifies  to 
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and  for  the  attracted  species  simplifies  to 
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The  average  yield  requires  a  numeric  integral. 
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For  the  attracted  species  use  the  substitution  u  =  exp 
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For  the  repelled  species  use  the  substitution  u  =  exp 
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The  numerator  is  calculated  numerically  using  200  points  and  Simpson’s  integration. 

3.1.2  Fontheim  Currents 

The  Fontheim  distribution  includes  both  a  Gaussian  and  a  power  law  term.  This  distribution  is 
only  applicable  for  electrons. 

Gaussian  Term 

The  Gaussian  integrals  are  slightly  more  complex  than  the  Maxwellian. 
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Use  the  substitution  u 
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where  L  is  now  max 


f-E0  —  4>  —  E0  ^ 

A  ’  A 

y  gauss  gauss  y 

The  substitution  x  =  exp(-  u2 )  can  be  used  to  solve  the  remaining  integral. 

j  =  ^qAgauss  yexp(-  L2)+  <q(4  +  E0)Agauss  ^erfc(L)  (82) 

The  average  yield  requires  a  numeric  integral.  Using  the  same  substitution  and  definition  of  L  we 
have 
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The  integral  in  the  numerator  can  be  divided  into  three  in  order  to  isolate  the  peak  at  u  =  0.  If 
L  >  0,  then  the  first  two  terms  are  zero  and  £,  =  L. 
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The  first  and  third  terms  can  be  evaluated  using  the  substitution  x  =  exp(-  u2 ).  Then 

u  =  ±V-  lnx  ,  where  the  sign  is  used  in  the  first  term  and  the  “+”  sign  in  the  third.  Under  the 
assumption  that  £,  is  small,  symmetry  can  be  used  to  write  the  second  tenn  in  terms  of  the  error 
function.  The  numerator  can  be  written  as 
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Power  Law 
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Using  the  substitutions  x  =  E-<j)  and  then  z  =x'(“'1)  and  z  =x'“  we  get 
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The  same  technique  is  used  to  evaluate  the  average  yield. 
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3.1.3  Kappa  Distribution  Currents 
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which  for  the  repelled  species  simplifies  to 
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(90) 


for  k  >  1 . 

For  the  attracted  species  the  integral  requires  a  substitution. 
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The  average  yield  requires  a  numeric  integral. 
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Use  the  substitution  u  =  1  + 


E  ±  cj) 
K0 


Jy(k0(u  -l)  +  (J))(k0(u  -l)  +  (|))(u)  k  1  du 

(Y)  =  ^ - - - 

|(K0(u-l)  +  ^)(u)“K_1du 
V 


(93) 


where  L'  =  1  ± 


K0 


for  the  repelled  species  and 


V 


for  the  attracted  species. 


Then  use  the  substitutions  y  =  uK_1  and  y=  uK  for  an  even  distribution  of  points. 


where  Lj  = 


1  +  2U 
K0 


^-D 


V 


and  L2  = 


1±A 

K0 


-1/ 


for  the  repelled  species  and  Li  =  L2  =  1  for 


the  attracted  species. 
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3.1.4  Convected  Maxwellian  Currents 

As  the  convected  Maxwellian  distribution  depends  on  the  velocity  of  the  plasma  with  respect  to 
the  spacecraft  and  the  angle  between  the  bulk  plasma  velocity  and  the  particle  velocity,  the  full 
angular  integral  is  needed. 


js  =  j,  (n,  O,  U)  =  q  JJJd3  v(  v  ■  h)fs  (v) 
The  coordinate  system  is  shown  in  Figure  8. 


(95) 


Z 


Figure  8.  Coordinate  System  Used  for  Convected  Maxwellian  Environment,  in  Spacecraft  Frame. 

Written  as  an  integral  over  the  velocity  at  infinity,  v®,  we  have 
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sin  cos  (96) 


S  J 


/2e<() 


where  the  lower  velocity  limit  L=0  for  an  attractive  potential  and  L  =  — 1  otherwise.  %  is  the 

V  ms 

angle  between  the  flow  vector  U  and  the  incident  velocity  at  infinity  v®.  %  is  dynamically 
mapped  to  the  incident  angle  S  and  is  given  by 


cos  x  =  cos  S^  cos  Su  +  sin  Sm  sin  Su  cos  cp . 


(97) 


where  Su  is  the  polar  angle  of  the  flow  velocity  vector  and  S®  is  the  polar  angle  of  v®.  With  this 
substitution,  the  azimuthal  integrand  has  the  form"  exp(-a  cos(cp)),  and  noting  that  the  zeroth- 
order  modified  Bessel  function  of  the  first  kind  Io  is  given  by 

2ftl0(a)=  J  d(pe~acosip ,  (98) 

0 


we  can  write  the  current  as 


*  The  evenness  of  the  (even-order)  Bessel  function  swallows  a  minus  sign  here. 
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The  velocity  and  polar  angle  integrals  are  performed  numerically.  The  angle  and  magnitude  of 
the  flow  velocity  vector  (defined  in  the  spacecraft  frame)  are  specified  by  the  user,  and  the  angle 
at  infinity  0ro  is  obtained  from  orbit  relations  in  a  1/r  potential,  as  a  function  of  the  total  energy 
mv2 

E  = - —  and  launch  angle  9  L.  From  reference  32  we  obtain^ 

2e 


<t> 


r  2(E  +  c|))sin2  9 


1  +  /l+4E(E2+(^sin29  cos(9l  -  90 ) 


<t> 


(100) 


Measuring  orbital  angles  from  the  polar  axis  that  is  normal  to  the  sphere  at  the  launch  point  and 
setting  a  =  r  determines  the  angle  90 
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Setting  r  =  oo(9  =  9oo)  gives 


008(9^  -  90)  =  -1 


4E(E  +  (j>)  2 


sin"  9 


(102) 


In  the  limit  as  U— >0  we  recover  the  expected  isotropic/uniform  flux  density,  while  as  the  Mach 
number  increases  the  flux  approaches  a  cosine  dependence  on  relative  particle  angle.33 

For  full  consistency,  the  secondary  electron  emission  current  due  to  ion  or  electron  impact  and 
the  backscattered  current  due  to  electron  impact  would  be  given  by 


55  Starting  from  Goldstein  eqn  3-46,  substitute  angular  momentum  l  =  rxp  =  rmv  sin(  9 ),  potential  energy  ®  =  k/r, 
and  l/2mv2  =  E  +  <t>,  then  for  a  launch  radius  r  =  a). 
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As  this  is  a  time  consuming  calculation  and  electron  Mach  numbers  are  always  under  1 ,  we 
compute  the  yields  due  to  incident  electrons  using  an  isotropic  Maxwellian  distribution. 

•sec, back  /-y  \  sec, back  .incident 

J  e  \  maxwellian  /  e  J  e  (104) 

The  secondary  emission  due  to  ions  is  computed  using  the  full  angular  integral. 

The  effective  density  of  each  component  is  the  expectation  value  of  the  distribution  function.  The 
techniques  used  to  compute  the  average  values  of  the  square  root  of  the  inverse  of  the  energy  are 
the  same  as  used  to  evaluate  the  current  integrals  above. 

3.1.5  Measured  Currents 


If  F(E)dE  is  the  measured  flux  of  ions  or  electrons  striking  a  zero  potential  surface  with  energy 
between  E  and  E+dE,  after  correction  for  the  actual  surface  potential  of  the  detector,  then 


j  =  q 


E  +  f 


F(E  ±  (())dE  =  q  J  [  ^-^Vujdu 


L±(j> 


U  ) 


(105) 


where  the  substitution  u=E±(|>  was  made. 

If  the  measurement  is  expressed  as  a  set  of  energy  bins  in  ascending  order  with  flux 
Fi(uj+i  -  uj )  in  bin  i  which  extends  from  Uj  to  U  j+| ,  the  integral  is  discretized  as 


r  fu  +  <|0 
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7u,ui+1  +4)A 
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where  uJ+i  is  the  lower  boundary  of  the  first  bin  entirely  within  the  domain  of  integration, 

Uj  <±<|)<Uj  +  1. 


The  average  yield  is  then  given  by 
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3.2  Space  charge  limiting  of  photoemission 

For  a  barrier  height  of  Ei,  the  escaping  photocurrent  for  a  surface  exposed  to  the  sun  at  an  angle 

MFun  is  given  by  the  formula  jpShC0^e(El )  =  'sunYsun  (Fi  )max(cos(\|/sun  ),0),  where  IsunYsun (E, ) 

is  the  yield  of  photoelectrons  with  energy  greater  than  Ei  from  the  surface  exposed  to  a  normally 
incident  solar  spectrum.  In  the  absence  of  electrostatic  saddle-point  or  space  charge  barrier 
formation,  the  barrier  to  escape  of  photoelectrons  is  the  surface  potential  for  positive  potentials 
and  zero  for  negative  potentials.  Photoelectron  space  charge  provides  an  additional  possibility  for 
increasing  the  barrier  in  excess  of  the  surface  potential  or  imposing  a  barrier  in  front  of  a 
negative  potential  surface.  If  the  user  requests  that  space  charge  limiting  of  the  photocurrent  be 
included  in  the  charging  calculation,  the  following  algorithm  is  used  to  determine  the  barrier 
height,  Ei. 

First,  if  the  ratio  of  the  available  photocurrent  to  the  plasma  density  is  less  than  2%  of  the  typical 

value  for  this  ratio  at  Earth  orbit,  ^photo  — IX  <  0.02  there  is  no  barrier. 

78xl(T6  na 


Otherwise,  a  “maximum  barrier,”  Bmax,  is  computed.  Bmax  is  the  lesser  of  the  barrier  calculated  in 
two  different  ways.  The  first  estimate,  B^ax ,  is  given  by  the  condition  that  the  barrier  occurs  no 
further  from  the  surface  than  the  distance  at  which  the  current  drops  by  a  factor  of  two,  i.e., 
r  f  .  <  2Rp  ,  where  the  symbols  represent  the  distance  of  the  barrier  from  the  local  center  of 

barrier  c  ’  J  r 

curvature  and  the  effective  radius  of  curvature  of  the  surface.  The  second  estimate,  B^,jx  ,  is 

given  by  the  condition  that  the  photoelectron  density  beyond  the  barrier  must  be  no  less  than  half 
the  ambient  plasma  density.  If  either  condition  gives  a  negative  value,  the  barrier  is  zero. 

13  max  is  found  using  an  analytic  fit  to  numeric  calculations  that  assume  only  that  the  shape  of  the 
photoemission  spectrum  is  the  default  shape  specified  in  Section  2.1.3.  Without  loss  of 
generality,  the  minimum  potential  is  taken  to  be  zero,  therefore  the  surface  potential  is  Ei.  In  a 
locally  spherical  geometry,  the  space  charge  between  the  surface  and  the  barrier  is  given  by 
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(109) 
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where  the  Heaviside  step  function  has  been  used  to  reflect  the  fact  that  electrons  with  energy  less 
than  Ei  both  leave  and  return,  while  those  with  energy  greater  than  Ei  only  leave.  The  surface 
potential  can  be  found  by  integrating  Poisson’s  equation  inward  from  the  barrier  to  the  surface. 
For  any  combination  of  escaping  photo  current  and  effective  radius  of  curvature  there  is  a  unique 
value  of  escaping  energy  Ei  for  which  the  surface  potential  is  also  Ei.  The  fit  to  these  self- 
consistent  barrier  calculations  used  in  Nascap-2k  is 

B<1  =  (24.772  +  13.755Rc  -7.6419R;)(j“r(0)f42’MI54"‘  (110) 


Bmax  is  found  using  an  analytic  fit  to  the  solution  of  equating  the  charge  density  of  the  escaping 
photocurrent  to  half  the  ambient  plasma  density.  Beyond  the  barrier  in  the  ambient  plasma,  the 
charge  density  of  the  escaping  photocurrent  is 


/»oo 


Equating  Pescaping  to  half  the  plasma  density  and  rearranging  gives 


/»00 


1  ena  |2e 

2j;r(o)  iK 


(in) 


(112) 


where  the  left  hand  side  depends  only  on  the  shape  of  the  photoemission  spectrum  and  the 
surface  potential,  and  the  right  hand  side  depends  only  on  the  problem  parameters  and  incident 
angle.  Again,  when  solving  this  equation,  the  shape  of  the  photoemission  spectrum  is  taken  to  be 
default  shape  specified  in  Section  2.1.3.  The  maximum  barrier  is  given  by  zero  surface  potential. 
For  zero  surface  potential,  the  solution  is  well- fit  for  barriers  up  to  12  volts  by 

•escape 

1 .436 x  1(T8  Jphot°  V  ’  =  1 0“3 (9.7 1  +  B® (0. 1 99  +  B®  (2.04  +  0.378 B^ )))  (11 3) 

c  n 


where  the  scaling  is  such  that  the  left  hand  side  is  unity  for  typical  1  AU  parameters.  To  simplify 
the  solution  of  the  cubic  equation  above,  we  obtain  an  initial  guess  by  first  solving  the  quadratic 
fit 


1.436x10  8 


en„ 


10-J(61.3-B™(39.7-9.41Bm)) 
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and  then  solve  the  cubic  using  Newton  iteration. 


Once  the  maximum  barrier,  B„viv  =  minlB'!'  ,  B':'  ),  has  been  detennined,  the  actual  barrier  Ei 
can  be  determined.  If  (j)s  >  Bmax  >  0 ,  the  space  charge  is  not  sufficient  to  cause  a  negative 
potential  in  front  of  the  surface,  so  the  barrier  is  equal  to  the  surface  potential. 

If  the  surface  potential  is  less  than  the  maximum  barrier  (Bmax  >  (j)s,  including  cases  with 
negative  surface  potential)  a  final  condition  is  used  to  detennine  the  barrier.  The  escaping 
photocurrent  must  be  sufficient  to  give  photoelectron  density  comparable  to  the  ambient  density 
beyond  the  barrier  in  the  ambient  plasma.  This  condition  is  approximated  as 

j“:r(E1)inieipE35  (115) 


If  Ei  =  Bmax  satisfies  this  inequality,  the  barrier  equals  the  maximum  barrier.  If  not,  the  equation 
is  numerically  solved  for  a  lower  barrier  Ei  that  gives  the  equality. 

3.3  Shadowing 

A  surface  element  is  shadowed  from  the  sun  if,  either 


(1)  Its  nonnal  points  anti-sunward,  or 

(2)  Its  centroid  lies  behind  the  projection  of  another  sun- facing  surface. 


The  detennination  of  the  second  is  handled  in  a  coordinate  system  in  which  the  z-axis  is  the 
direction  to  the  sun.  The  unit  vectors  along  the  axes  of  the  rotated  coordinate  system  are  given  by 
sx,  sy,  sz.  The  position  of  a  surface  element  along  the  z-axis  in  the  rotated  system  is  given  by 


rz  =  —  ^(sz  •ri),  where  r,  are  the  nodes  of  element  and  N  is  the  number  of  nodes.  The 
N  j 

projection  of  a  surface  element  on  the  x-y  plane  of  the  rotated  system  is  given  by  (s  x  •  r,  s  y  •  r  ) . 


Surface  A  is  shadowed  by  surface  B  if,  rz(A)  <  rz(B)  and  the  centroid  of  the  projection  of  surface 
A  is  inside  the  projection  of  surface  B. 

3.4  Implicit  Charging  Using  the  Boundary  Element  Method 

Among  the  difficulties  of  developing  accurate  and  robust  algorithms  for  spacecraft  charging  has 
been  the  inability  to  calculate  electric  fields  accurately,  let  alone  to  predict  how  electric  fields 
will  change  as  a  result  of  surface  potential  changes.  Nascap-2k  uses  the  Boundary  Element 
Method  (BEM)34  to  calculate  accurate  electric  fields  and  as  the  basis  for  implicit  charging 
equations. 

3.4.1  Boundary  Element  Method  Algorithm 

The  Boundary  Element  Method  is  a  means  for  relating  fields  and  potentials  in  a  region  to  sources 
on  the  boundary  of  the  region.  It  is  comparable  to  a  sum  over  the  coulomb  field  of  all  the  charges 
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in  a  region  rather  than  an  iterative  field  solution.  In  our  case,  the  region  is  the  space  exterior  to  a 
spacecraft,  and  the  boundary  is  the  spacecraft  surface.  Also,  we  assume  the  “free  space  Green’s 
function”,  i.e.,  the  potentials  in  the  region  obey  Laplace’s  equation. 

The  sources  are  sheets  of  charge  coincident  with  the  spacecraft  model’s  surface  elements.  We 
assume  that  each  surface  element,  j,  has  a  constant  charge  density,  a7.  The  familiar  relation  for 
the  potential  of  a  point  charge  then  generalizes  to  an  integral  over  the  object  surface: 


<t> 


q 

47tsor 


(116) 


where  (j);,  which  could  be  the  potential  at  any  point  in  space,  is  considered  the  potential  at  the 
center  of  a  surface  element.  Similarly,  the  familiar  relation  for  the  electric  field  of  a  point  charge 
generalizes  to: 


E  =  ■ 


4ti£  r 


-A 


(47t£0)El=^[d!rA(r,-rJ) 

j  J  hj 
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where  E;  is  the  electric  field  at  some  point  in  space,  the  point  again  taken  at  the  center  of  a 
surface  element. 


We  express  the  relations  of  potential  and  electric  field  to  charge  density  as  matrices: 

♦ific-'K  'Km  (118) 

These  can  be  combined  to  obtain  a  relation  between  nonnal  electric  field  and  voltage: 

(E  -  n -FikCkjVJ  (119) 

This  last  relation  is  the  key  to  developing  relations  between  surface  charge,  surface  potential,  and 
surface  currents  in  order  to  derive  charging  equations. 

3.4.2  Doing  the  Integrals 

To  get  C'1  we  need  to  do  integrals  of  the  form 


r 


J 


(120) 


There  are  tricks  to  doing  these  integrals,  which  can  be  found  in  the  literature.  Denote  the  field 
point,  r;  as  P,  and  take  the  domain  of  integration  as  triangle  ABC.  Then,  the  vector  from  the  field 
point  to  anywhere  on  the  triangle  can  be  parameterized  as 

r;j  =  PA  +  uAB  +  uvBC  (121) 
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where  PA,  AB,  and  BC  are  vectors  between  pairs  of  points,  and  u  and  v  are  parameters,  each  in 
the  interval  [0,1].  The  square  of  the  distance,  pj,  then  becomes 

r^rrriJ=  ZTk/uV  (122) 

0<k./<2 

where  the  coefficients  Ty  are  formed  from  pairwise  scalar  products  of  the  three  vectors  above. 
Now,  the  integral  can  be  expressed  as 
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where  the 
The  inner 


Vi  coefficients  are  functions  of  v. 

integral  (over  u)  may  be  found  in  standard  integral  tables 
1  Vc  +  bx  +  ax2 
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We  are  left  to  do  the  outer  integral  (over  v)  numerically.  To  facilitate  this,  we  select  the  vertex  A 
such  that  the  scalar  product  PA-BC  is  the  minimum  of  the  three  choices.  Then,  very  few 
integration  points  are  needed  in  the  outer  integral  (over  v).  (We  use  a  5-point  Simpson’s  rule  in 
the  present  implementation.) 

The  integrals  for  the  electric  field  are  similar,  although  there  are  more  of  them.  The  same 
techniques  apply. 

3.4.3  Test  for  Accuracy 

We  tested  the  algorithm  for  accuracy  by  calculating  the  electric  fields  on  the  surface  elements  of 
a  unifonn  potential  sphere.  The  sphere  model  in  Figure  9  was  originally  defined  using  Patran, 
and  has  90  surface  elements  and  92  nodes.  It  provides  a  fairly  coarse  mesh  over  the  surface  of 
the  sphere. 
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Figure  9.  Sphere  Defined  Using  Patran. 

The  next  figure  shows  the  result  of  the  BEM  calculation.  Note  that  the  total  variation  of  electric 
field  over  the  surface  is  only  three  percent.  Also,  the  apparent  radius  of  the  sphere  (given  by  V/E) 
is  within  three  percent  of  the  nominal  radius. 
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Figure  10.  Normal  Component  of  Electric  Field  at  Surface  of  Sphere  as  Computed  Using  the  BEM 
Approach. 
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By  contrast,  the  same  calculation  using  Nascap-2k ’s  finite  element  algorithm  with  higher  order 
elements  gives  an  electric  field  variation  of  seven  percent,  and  a  value  of  electric  field  that  is 
about  fifteen  percent  high.  A  finite  element  algorithm  with  linear  elements  gives  an  eight  percent 
variation,  and  a  value  that  is  about  twenty-five  percent  high. 

3.4.4  Implicit  Charging 

To  develop  charging  equations  we  need  to  express  physical  charges  on  physical  surfaces  in  terms 
of  the  voltages  on  the  same  objects.  We  can  then  proceed  to  calculate  what  voltage  changes  will 
be  produced  by  changes  in  charge  (currents).  Because  the  interior  of  a  spacecraft  is  not  free 
space,  the  physical  charge  densities  bear  no  relation  to  the  charge  density  a’s  used  in  the 
derivation  of  the  relation  between  the  normal  electric  field  and  the  voltage  of  Boundary  Element 
Method. 

To  get  the  physical  charges  we  use  Maxwell’s  divergence  equation  in  the  form  a  =  V-D.  Figure 
1 1  shows  the  “Gaussian  pillboxes”  used  to  calculate  the  actual  surface  charges  on  insulating 
surfaces  and  on  conductors,  as  described  below. 


Figure  11.  “Gaussian  Pillboxes”  Used  to  Calculate  the  Actual  Surface  Charges  on  Insulating 
Surfaces  and  on  Conductors. 


For  an  insulating  surface  (upper  pillbox  in  Figure  1 1),  the  external  field  is  E-n,  which  we  obtain 
from  the  BEM  solution.  The  internal  field  is  related  to  the  capacitance  between  the  insulating 
surface  and  its  underlying  conductor.  Thus,  the  total  charge,  qt  on  such  a  surface  is  given  by: 


q,  =  A^En),  +Ci0(<t>, -<I>J 


(125) 


where  A  is  the  surface  area.  For  a  conductor,  since  charges  are  mobile,  it  is  not  useful  to  know 
the  charge  on  each  individual  surface,  but  we  need  to  work  with  the  total  charge  on  the 
conductor.  Conducting  surfaces  include  the  obvious  “bare”  surfaces  (lower  pillbox  in  Figure  11), 
as  well  as  the  surfaces  that  underlie  the  insulating  surfaces  (middle  pillbox  in  Figure  1 1).  In  both 
cases,  we  have  zero  electric  field  internal  to  the  metal.  To  obtain  the  total  charge  on  all  the 
surfaces  of  the  conductor,  we  need  to  sum  the  external-field  charge  terms  for  the  bare  conducting 
surfaces,  plus  the  capacitive-charge  terms  for  the  insulator-metal  interfaces: 

Q,=£a,(E-ii),-  Ecjfc-G 

bare  insulators  ( 1  76) 
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We  previously  found  the  BEM  expression  for  the  external  fields  in  terms  of  the  cell  potentials: 

(JL  -  n  )i  =  FikCkji]>j  (127) 

Substituting  this  into  the  equations  for  q,  and  Qc,  perfonning  the  indicated  sums,  and  adding  the 
capacitive  terms,  we  get  the  matrix  equation: 


Q=GO>  (128) 

where  the  vectors  are  composed  of  contributions  from  all  the  insulating  surfaces  and  a  tenn  for 
each  conductor  representing  all  the  bare  conducting  surfaces  of  that  conductor.  Conductors  that 
are  biased  with  respect  to  another  conductor  are  treated  as  a  single  conductor  at  this  point  in  the 
calculation. 


Q=  {qi,q2,  ...  ,qn,  Qc}  (129) 

®  =  {<t>l,  Cj>2,  •••,  <j>n,  (130) 


where  Q  and  O,  the  charges  and  potentials,  are  related  by  the  charging  matrix,  G. 

We  are  looking  for  an  equation  that  relates  currents  to  voltage  changes.  So,  we  differentiate  the 
charge  equation  in  time: 


i=q=go 


(131) 


Discretize  to  a  finite  time  interval: 


IAt  =  G[<J>(t  +  At)-<J>(t)]  (132) 

Implicitize  by  evaluating  current  at  the  final  time  (also,  for  simplicity,  changing  [t+At,t]  to  [t,0]): 

l(t)t  =  G[a>(t)-O(0)]  (133) 

Linearize  currents  with  respect  to  voltage  (since  we  do  not  know  the  final  voltages  at  which  the 
current  is  to  be  evaluated): 


i(t)»i(o)+r[®(t)-®(o)] 


(134) 


And,  solve  the  equation: 


[G  -  I't]o(t)  =  [G  -  I't]<J>(0)  +  l(0)t  (135) 

This  gives  us  a  straightforward  matrix  equation.  Everything  on  the  right-hand-side  is  known,  and 
we  can  solve  using  standard  linear  algebra  equation  solver  packages. 

Any  user  specified  values  for  the  bias  of  one  conductor  with  respect  to  another  is  reapplied  each 
timestep.  The  v  x  B  contribution  to  the  surface  potential  is  also  reapplied  each  timestep. 
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Before  proceeding  to  examples,  it  is  worth  commenting  on  the  derivative  of  current, 

Is  =  SI  ,/a*,- 

Consider  three  cases: 


1 .  For  current  sources  such  as  incident  plasma  current,  we  usually  approximate  the  current 
as  a  function  of  the  local  voltage  only.  This  gives  a  diagonal  term  in  II .  Since  such 

currents  decay  exponentially,  some  care  must  be  taken  not  to  underestimate  the  change  in 
current  if  the  voltage  is  changing  in  such  a  direction  that  a  component  of  current  is 
increasing  (i.e.,  electron  current  for  a  surface  whose  potential  is  increasing  (toward  zero) 
from  a  large  negative  potential). 


For  example,  suppose  we  have  (j)  =  I((j))  =  1  -  and  c))(0)=-2.  The  implicit  solution  for 
t()(t)  is  4>(t)  =  ^(0X1  1 0^+ |p  we  use  tjie  aciuai  derivative  I'  =  -exp((|)(o)),  we 

obtain  the  upper  curve  in  Figure  12,  which  works  satisfactorily  for  sensible  timesteps  like 
0.1  or  1,  but  overshoots  the  solution  ((j)(oo)  =  0)  badly  for  long  timesteps  like  10  or  1000. 
However,  if  we  recognize  that  zero  is  a  value  we  wish  not  to  overshoot  and  accordingly 

set  I'  =  ^ — XfoCO))  we  013^^  the  lower  curve  in  Figure  12,  which  gives  equally  good 

0-(j)(0) 

results  for  short  timesteps  and  does  not  overshoot  the  solution,  even  for  a  timestep  of 
1000. 


0.1  1  10  100  1000 

Timestep 


Figure  12.  Example  Problem  Solutions  for  Short  and  Long  Timesteps  for  Two  Choices  of  Current 
Derivative. 


2.  Conductivity  current,  such  as  I;  =  cr(c()c-ct)i),  contributes  off-diagonal  as  well  as  diagonal 

ar 


terms  to  the  current  derivative  matrix, 
contributes  many  more  off-diagonal  terms. 


ai, 

va<t>i 


-ct; 


34k 


■  +CT 


.  Surface  conductivity 
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3.  The  case  that  is  particularly  problematic  is  I;  =  T(E-n).  This  occurs  for  a  photo-emitting 
surface  at  negative  potentials.  The  problem  is  that  the  local  electric  field  is  a  function  of 
the  potentials  on  all  the  surfaces.  However,  the  BEM  provides  us  with  exactly  that 
function.  We  can  now  compute  the  term 

ly  =  SI; /d(t>j  =  SIj  / SEj  x  <3Ej / Scjjj  (136) 

using  the  relation 

(E  -  n),  =FikCljOj  (137) 


derived  from  the  BEM. 

In  Nascap-2k  we  impose  some  additional  conditions  on  the  current  derivative  in  order  to  promote 
stability,  some  of  which  are  related  to  the  sharp  cutoffs  that  occur  with  secondary  and 
photoemission  currents.  The  problem  is  illustrated  by  Figure  13.  If  we  charge  negative  from 
point  A  and  use  the  actual  derivative  we  fail  to  anticipate  the  sharp  increase  in  secondary 
emission  current  and  overshoot  negative.  If  we  charge  positive  from  point  B  and  use  the  actual 
derivative  we  fail  to  anticipate  the  sharp  drop  in  secondary  emission  current  and  overshoot 
positive.  If  we  charge  negative  from  point  C,  then  the  actual  derivative  predicts  unphysically 
enhanced  secondary  emission  current  at  negative  potentials,  with  the  result  that  negative 
charging  becomes  glacially  slow.  The  rules  summarized  below  are  designed  to  promote 
reasonable  behavior  for  charging  of  surfaces  encountering  suppressed  emission  current. 


Figure  13.  Notional  Curve  of  Effective  Secondary  Emission  Current  versus  Surface  Potential. 

1 .  If  the  secondary  electrons  are  limited  by  positive  potential  or  electric  field  barrier  and  the 
net  current  would  be  positive  were  they  not  limited,  then  the  rate  of  increase  of  emission 
current  with  drop  in  potential  is  at  least  doubled,  and  for  potential  above  2  V  the  rate  is 
adjusted  such  that  emission  is  turned  on  fully  should  the  potential  drop  to  zero.  This 
prevents  (or,  at  least,  minimizes)  oscillations  between  high  positive  and  high  negative 
potential  when  the  true  current  balance  occurs  at  near  zero  potential. 
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2.  For  the  Fontheim  distribution,  if  the  electric  field  is  positive  and  the  surface  has  negative 
current  before  any  limiting  of  low  energy  currents,  then  the  derivative  is  set  to  ignore  the 
increase  in  emission  due  to  decreasing  potential.  This  allows  the  surface  potential  to 
decrease  at  a  reasonable  rate.  Otherwise  the  low  energy  currents  would  go  from 
suppressed  to  unphysically  enhanced  and  prevent  the  correct  negative  charging. 

3.  If  the  electric  field  is  negative  and  the  surface  has  positive  current  before  any  limiting  of 
low  energy  currents,  then  the  derivative  is  set  to  turn  off  the  emission  when  the  potential 
reaches  2  V.  Otherwise,  the  suppression  of  emission  at  positive  potential  would  not  be 
anticipated  and  the  unsuppressed  emission  would  drive  the  surface  far  above  its  true 
floating  potential. 

4.  If  a  surface  is  negative  but  has  positive  current  and  positive  nonnal  electric  field,  the 
derivative  is  set  so  that  the  current  is  turned  off  after  rising  to  a  potential  of  no  more  than 
+5  V.  This  prevents  charging  from  negative  to  unphysical  positive  values  when  emission 
is  already  suppressed. 

5.  The  current  derivative  is  set  to  allow  a  potential  change  of  no  more  than  2  kV. 

6.  The  current  derivative  is  set  to  allow  a  potential  change  of  at  least  the  minimum  of  1  V 
and  one-one  hundredth  of  the  electron  temperature. 

3.4.5  Applicability  at  Finite  Debye  Length 

The  effect  of  using  the  BEM  capacitances  to  calculate  charging  dynamics  is  that  the  “to  space” 
capacitance  is  underestimated,  and  thus  the  “overall”  (to  space)  charging  rate  of  the  spacecraft  is 
too  fast.  For  short  Debye  length  and  low  temperature  the  capacitance  is  inversely  proportional  to 
the  Debye  length,  whereas  the  BEM  (vacuum)  capacitance  is  inversely  proportional  to  the 
spacecraft  size.  For  typical  cases  the  charging  rate  may  be  too  fast  by  a  factor  of -100.  Note  that 
even  the  “correct”  charging  rate  is  fast  compared  with  the  differential  charging  rate. 

Two  factors  that  mitigate  the  charging  rate  issue  outlined  above  are; 

1 .  At  high  potentials  the  capacitance  is  inversely  proportional  to  the  sheath  size  rather  than 
to  the  Debye  length,  so  that  typical  errors  in  charging  rate  are  factors  of  ~2  to  10  rather 
than  —100; 

2.  Typically,  in  Nascap-2k  calculations  the  timestep  is  far  too  long  to  resolve  the  vacuum 
charging  rate,  so  that  the  charging  rate  in  the  simulation  is  detennined  by  the  choice  of 
timestep  rather  than  the  actual  capacitance. 

Note  that  the  rate  of  differential  charging  is  not  affected  by  the  plasma,  as  the  capacitances 
governing  differential  changing  depend  on  material  thicknesses  and  electronic  component  values 
rather  than  on  spacecraft  size  and  plasma  effects.  Also,  the  quasi-steady-state  and  equilibrium 
potentials  are  not  affected,  because  these  depend  on  currents  which  in  turn  depend  on  potentials 
and  electric  fields.  The  electric  fields  used  to  calculate  currents  are  calculated  from  the  finite 
element  potential  solution,  and  therefore  do  include  the  effects  of  a  finite  Debye  length  plasma. 

3.4.6  Example:  Sunlit  Sphere 

To  illustrate  the  implicit  charging  model,  we  calculate  the  charging  of  a  sunlit  Teflon  sphere  in  a 
1  cm'3,  20  keV  plasma.  The  NASCAP/GEO  version  of  this  was  published  in  1978. 36  A  version  of 
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the  original  result  is  shown  below  in  Figure  14.  The  sun  direction  is  (1,1,1)  (from  the  upper 
right).  The  dark  side  of  the  sphere  gradually  charges  negative  due  to  incident  plasma  electrons, 
while  photoemission  grounds  the  sunlit  side.  Eventually,  however,  the  strong  negative  potentials 
due  to  the  dark  surfaces  wrap  around  the  sphere  and  form  a  barrier  to  photoelectron  escape.  The 
potential  of  each  sunlit  surface  is  subsequently  detennined  by  the  condition  that  its  electric  field 
has  a  small,  positive  value  (seen  in  Figure  14  by  the  coarse  contour  spacing  near  the  upper  right 
spacecraft  surface),  so  that  just  the  right  fraction  of  photoelectrons  can  escape  over  the  barrier. 
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Figure  14.  Potentials  About  Sunlit  QSphere  in  Space  as  Computed  By  NASCAP/GEO. 

Figure  15  and  Figure  16  show  the  BEM  solution  of  the  same  problem  for  a  PATRAN  sphere  and 
a  QSphere  respectively.  The  point  on  the  sphere  furthest  from  the  sun  is  the  least  negative  point. 
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Figure  15.  Potentials  on  Sunlit  PATRAN  Sphere  in  Space  Viewed  From  Direction  (1,2,3). 
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Figure  16.  Potentials  on  Sunlit  QSphere  in  Space  Viewed  From  Direction  (1,2,3). 

Figure  17  is  another  view  of  the  BEM  solution,  showing  more  of  the  dark  side.  Despite  the  fairly 
coarse  gridding  on  the  sphere,  the  gradual  potential  gradient  on  the  sunlit  side  and  the  constant, 
large  negative  potential  on  the  dark  side  are  clearly  seen  in  the  BEM  solution. 
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Figure  17.  Potentials  on  Sunlit  PATRAN  Sphere  in  Space. 

3.5  Finite  Element  Method  to  Compute  Potentials  in  Space 

In  calculating  the  potential  in  three  dimensions  around  an  arbitrary  object,  Nascap-2k  employs  a 
finite  element  approach  using  right  parallelepiped  elements  and  nonlinear  edge  interpolants. 
Arbitrarily  nested  subdivision  (up  to  seven  levels)  is  used  to  resolve  important  object  features 
while  including  a  large  amount  of  space  around  the  spacecraft  for  determining  accurate  particle 
orbits  and  keeping  memory  requirements  reasonable.  Figure  18  shows  an  example  of  a 
computational  grid  with  progressively  finer  subdivision  to  give  fine  resolution  near  an 
experiment. 
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Figure  18.  Computational  Grid  With  Subdivision. 


Nascap-2k  solves  Poisson’s  equation 


-  V2(j)  =  p/s0 


by  solving  the  associated  variational  principle 


o 

II 

j  |  <yi 

r 

dV 

w  IJ 

L2  £o  J 

+  J  (j)V(j)  •  dS 


(138) 


The  first  term  in  the  volume  integrand  corresponds  to  the  Laplacian  operator;  the  second  term  is 
the  space  charge  contribution.  The  surface  contribution,  given  by  the  last  term,  is  used  to  obtain 
the  surface  electric  field. 

Poisson’s  equation  is  solved  with  surface  potentials  fixed.  To  determine  surface  field  values  for 
fixed  potential  elements,  we  use  the  fact  that  the  surface  element  electric  field  is  conjugate  to  its 
potential,  i.e.,  when  the  full  conjugate  gradient  matrix  multiplication  is  performed,  the  surface 
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element  electric  field  (corresponding  to  the  existing  potentials  and  fields)  is  the  residual 
corresponding  to  the  surface  element  potential.  This  is  an  improvement  over  capacitance  matrix 
methods  in  that  it  takes  full  account  of  the  nonlinear  space  charge  in  the  external  space.  In  the 
variational  calculation,  we  use  locally  defined  basis  sets,  that  is,  interpolants  within  each  cubic 
volume  element.  Fine  mesh  volumes  are  given  the  correct  variational  weight,  ensuring  the 
maintenance  of  accuracy  through  mesh  transition  regions.  The  potential  and  electric  field  are 
defined  at  each  grid  node.  The  potential  inside  each  volume  element  is  interpolated  from  the 
values  at  each  of  its  vertices. 


The  volume  integrals  are  solved  volume  element  by  volume  element.  The  potential  (and  electric 
field)  at  each  point  in  the  element  is  the  sum  of  interpolants  times  the  values  at  the  nodes. 


Tvtf  +  pfidv =  W  f  (vtfav,  +  y  f  ^dv. . 

^  So  ^  e  J  e  J  S0 


(139) 


The  potential  inside  each  element  is  given  by 

f(x,y,z)  =  ^Ni(x,y,z)(|)i 


(140) 


where  “i”  are  the  nodes  of  the  element  indexed  by  “e”,  and  the  N)  are  the  interpolants.  Note  that 
each  element  has  four  interpolants  per  node,  whose  coefficients  are  the  potential  and  three 
components  of  electric  field  at  the  node,  as  discussed  in  Section  3.5.1. 

We  then  have 


Vf  (x,y,z)  =  ^VNi(x,y,z)(j)1 


(141) 


which  allows  us  to  write 


WdV  =  W  J  XSVN11.(x>y,z)VNjp(x>y>z)t4,jpdV<  =  WXZW'«i4A 


ia  jp 


ia  jp 


(142) 


Where  the  indices  (a, (3)  run  over  the  potential  and  three  electric  field  components  at  the  nodes 
and  their  corresponding  interpolants.  Note  that  the  WjaJp  depend  only  on  the  shape  of  the  volume 
element. 


The  second  term  can  be  written 


P4> 


dv  =  X  f  dV,  =  y  f  y  (x.  y,  z)MV,  = 2>  W 


The  minimization  can  then  be  written  as  the  matrix  equation 


(143) 
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(144) 


8 

S(j) 


e  z  i  j  e  i 


8  i  Q  . 

-  -<pW(p  +  -<p  =0 

8<t>  l2  £o  J 


This  function  is  then  expressed  in  terms  of  the  vector,  <p,  of  all  nodal  values  and  a  sparse, 
symmetric,  positive  definite  matrix  W.  The  conjugate  gradient  solution  to  this  problem  requires 
that  we  be  able  to  form  the  matrix-vector  product  W<p,  which  we  can  do  by  adding  up  the 
contribution  to  W<p  from  each  finite  element  without  ever  forming  (or  storing)  the  sparse  but 
extremely  large  matrix  W. 

The  resulting  system  of  linear  equations  is  solved  using  the  Conjugate  Gradient  approach. 

3.5.1  Continuous  Field  Interpolants 

Nascap-2k  uses  a  finite  element  interpolation  scheme  that  has  strictly  continuous  electric  fields. 
This  approach  insures  the  continuity  of  the  electric  field  used  to  compute  particle  trajectories. 
Tri-quadratic  interpolants  were  also  considered,  but  experience  showed  that  the  difference  in 
weight  between  the  comer  nodes  and  the  edge-center  nodes  leads  to  a  poorly  conditioned  matrix. 

The  basic  interpolation  functions  consist  of  functions  (representing  nodal  potentials)  that  have 
unit  value  and  zero  slope  at  their  home  nodes,  and  zero  value  and  slope  at  opposite  nodes: 


(x-1)2 

1  -2x  +  2x2 


l-2x  +  2x2  (145) 


and  functions  (representing  nodal  electric  fields)  that  have  zero  value  and  unit  slope  at  their 
home  nodes  and  zero  value  and  slope  at  opposite  nodes: 


G0  =x(l-x)F0 
G,  =x(x-l)F, 


(146) 


The  coefficient  of  F  corresponds  to  the  potential  at  the  home  node,  and  the  coefficient  of  G 
corresponds  to  the  (negative  of  the)  electric  field  at  the  home  node.  (In  the  element  interior,  both 
F  and  G  are  needed  to  calculate  potential  and  potential  gradient.) 

The  functions  Fo  and  Go  are  shown  in  Figure  19.  It  can  be  verified  that  these  functions  can 
exactly  reproduce  constant,  linear,  and  quadratic  potentials. 
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Figure  19.  Interpolation  Functions  F0  and  G0. 

We  generalize  to  three  dimensions  by  assigning  to  each  comer  of  the  unit  cube  four  interpolation 
functions  corresponding  to  the  potential,  x-component  of  potential  gradient,  y-component  of 
potential  gradient,  and  z-component  of  potential  gradient  for  that  node.  The  interpolants 
corresponding  to  node  ijk  are  then  given  by 

NSk(x>y>z)  =  Fi(x)Fj(y)Fk(z) 

N|jk(x,y,z)  =  GI(x)FJ(y)Fk(z) 

(147 

Ni(x5y5z)  =  F,(x)GJ(y)Fk(z) 

N5k(x5y5z)=F,(x)Fj(y)Gk(z) 

The  coefficients  of  the  N  interpolants  are  the  node  potentials,  and  the  coefficients  of  the  N  ’  ’ 
interpolants  are  the  negative  Cartesian  components  of  the  nodal  electric  fields.  We  elect  to  omit 
the  higher  order  terms  of  the  form  GFG  or  GGG. 


The  potential  inside  each  element  is  given  by 

(Kx,y,z)  =  2N“k(x,y,z)(j)“k 

ijka 


and  the  components  of  the  electric  field  are  given  by 

E(x,  y,  z)  =  £  ^“kVN“k  (x,  y,  z) 

ijka 


(148) 


(149) 


3.5.2  Interpolating  Potentials  and  Fields  for  Special  Elements 

Special  elements  are  those  elements  which  contain  surfaces.  Because  surface-containing 
elements  are  not  empty  cubes,  the  potential  interpolation  functions  described  in  the  previous 
section  cannot  be  straightforwardly  applied.  To  develop  finite  element  matrices  for  these 
elements,  we  require  a  prescription  for  calculating  the  potential  and  electric  field  in  the  element 
interior.  Note  that  a  special  element  is  bounded  by  three  types  of  surfaces:  (1)  square  surfaces 
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bounded  by  grid  edges  and  shared  with  adjacent  (presumably  empty)  elements;  (2)  object 
surfaces;  and  (3)  surfaces  bounded  by  both  object  points  or  edges  and  grid  points  or  edges.  On 
type  (1)  surfaces  we  use  the  potential  and  field  interpolation  described  in  the  previous  section. 

On  object  surfaces,  the  potential  and  electric  field  must  be  expressed  in  terms  of  the  object’s 
surface  element  potentials  and  nonnal  fields.  Type  (3)  surfaces  must  smoothly  blend  between  the 
two.  We  describe  below  an  algorithm  to  express  the  potential  and  field  at  any  point  in  the 
element  volume  in  terms  of  the  grid  point  potentials  and  fields,  and  the  surface  element 
potentials  and  normal  fields.  The  algorithm  has  the  property  that,  when  applied  to  an  element 
with  six  type  (1)  surfaces,  it  reproduces  the  empty  cube  result  exactly. 

Let  R  be  a  point  in  a  volume  bounded  by  a  set  of  surfaces  { S } ,  each  of  which  may  be  a  triangle 
or  a  planar  convex  quadrilateral.  For  a  given  S,  let  P  be  the  point  on  S  nearest  R,  let  4>(P)  and 
E(P)  be  the  potential  and  electric  field  at  P,  and  n  be  the  unit  nonnal  (from  the  surface  point  P 
into  the  volume)  to  the  surface.  Let 


d  =  R  -  P 
d2  =  d  •  d 

Ks  =  As  /  d2  for  d  •  n  >  0  (150) 

N  =  'V'  Kg 

A) 

where  As  is  the  area  of  surface  S  and  Ks  is  its  weighting  function. 

Let  /  be  the  distance  from  P  in  direction  n  to  the  next  surface  intersection.  (In  practice,  it  is 
adequate  to  extend  /  to  the  intersection  with  the  surface  of  the  rectangular  parallelepiped 
bounding  the  volume.)  Then  the  contribution  of  S  to  the  potential  at  R  is 

(t>s  (R)  =  (Kg  /  N)[<KP)  -  (1  -  d  .  n  /l)(d  •  n)(E(P)  •  n )]  (151) 

and  the  total  potential  is  found  by  summing  the  contributions  from  all  surfaces. 

The  electric  field  is  found  by  differentiating  the  above.  The  contribution  of  S  to  the  electric  field 
is 


Eg  =  4>S[VN  /  N  -VKS  )  /  Ks  )]+  (Ks  /  N)[-  V^(P)  +  (1-  2d  •  n  /  l)n(E(P)  •  n)]  (152) 

Where  V(|)(P)  is  taken  tangential  to  the  surface,  and  is  calculated  using  the  continuous  field 
interpolants  if  S  is  a  full  square  in  empty  space,  and  otherwise  using  linear  interpolation  for 
triangles  and  bilinear  interpolation  for  quadrilaterals. 

This  approach  allows  us  to  interpolate  surface  fields  and  potentials  on  bounding  surfaces  in  a 
manner  that  maintains  the  continuous  electric  field  property.  Potentials  and  nonnal  field 
components  are  defined  at  the  centroids  of  the  original  object  surfaces  and  area- weighted 
averaged  at  surface  comers.  Potentials,  nonnal  fields,  and  tangential  field  components  at  all 
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surface  points  that  are  vertices  of  bounding  surfaces  are  expressed  as  linear  combinations  of 
centroid  potentials  and  fields. 


3.5.3  Boundary  Conditions 

Most  often  we  use  a  zero  potential  boundary  condition  at  the  outer  boundary  of  the 
computational  grid.  For  some  problems,  specifically,  Laplacian  potential  and  linearly  Debye 
screened  potentials,  a  different  approach  is  desirable.  For  these  problems,  we  apply  boundary 
conditions  by  extending  the  problem  beyond  the  computational  box  using  external  elements 
whose  inward  boundary  is  a  quarter-boundary-surface-element  (QBSE)  containing  one  problem 
node  and  which  extends  radially  outward  from  the  center  of  the  computational  box  (to  which  we 
henceforth  refer  as  “the  origin”).  (See  Figure  20.)  These  elements  are  added  to  the  sum  over  all 
volume  elements.  We  ignore  potential  variation  over  the  QBSE,  assuming  it  to  have  the  potential, 
^node,  of  its  problem  node,  so  that  this  element  makes  a  diagonal  contribution  to  W.  Within  that 
approximation  the  new  external  element  subtends  solid  angle  Q=  lAx2/R2,  where  the  QBSE  is 
nonnal  to  the  x-axis  at  distance  x  from  the  origin,  and  R  is  the  distance  of  the  center  of  the  QBSE 
from  the  origin.  (R  and  x  are  measured  in  outer  grid  units.)  The  potential  in  the  element  is 


<Kr)  =  ^  node 


'rQ 
V  r  J 


cxp(-(r-  R)//.n), 


(153) 


for  Debye  length  A,d-  We  can  now  proceed  to  evaluate  the  needed  integrals.  For  the  Laplacian 
(XD=oo)  case  we  get 


f  i(v*)2dv  =  i*L,.n 

/»00 

r2dr 

f-R" 

2 

J  2 

«/ 

R 

(154) 


so  that  after  properly  taking  account  of  the  mesh  spacing,  L,  the  external  element’s  contribution 
to  the  matrix  element  is 


Wvv  =}L(x/R2) 

where  x  and  R  are  measured  in  units  of  L. 


(155) 
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Figure  20.  An  Extended  Element  With  Analytically  Extrapolated  Potential  is  Appended  to  Each 
Boundary  Quarter-Square  to  Implement  Non-Zero  (Monopole)  Potential  Boundary  Conditions. 

For  the  case  of  finite  Debye  length,  the  evaluation  contains  more  terms  and  involves  integrals 
that  cannot  be  expressed  in  closed  form,  so  we  approximate 


dr 

d  R 


e 


-r/ X 


r 


n 


"R+A, 

dr 

J  R 


to  eventually  get 


W 


vv={L(x/R2)  ^  +  2pln 


1  +  — 
2p  j 


+  - 


1  ,  (L)3x 

1  +  2p  j  8pA2 


(156) 


(157) 


where  p =R  LA,. 


It  is  recommended  that  zero  potential  boundary  conditions  be  used  if  A  is  comparable  to  or 
shorter  than  an  outer  grid  mesh  unit,  and  for  most  particle  tracking  problems  for  which  no  clear 
value  exists  for  the  Debye  length. 


3.6  Space  Charge  Stabilized  Poisson  Iteration 

Poisson’s  equation  can  be  written  dimensionlessly  as 

-V2®  =  (n,-n,)/AJ 


(158) 


where 


0  =  -^,andA2=^^  =  s0  — ^  (159) 

0  L2  N0eL 

A  is  the  dimensionless  Debye  length,  N0  is  the  ambient  plasma  density,  n,  =  Ni/N0,  ne  =  Ne/N0, 
and  the  Laplacian  is  also  nonnalized  by  L2. 
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The  traditional  approach  to  the  solution  of  Equation  (158)  has  been  an  explicit  iteration  of  the 
fonn 


-V2®v  =A~2(ni(®v-1)-ne(®v-1))  (160) 

where  v  is  the  iteration  index  and  the  charge  density  is  determined  using  the  potentials  of  the 
previous  iteration.  This  method  can  be  shown  to  be  unstable  when  the  Debye  length,  A-d, 
becomes  small  with  respect  to  other  scale  lengths  of  the  problem.  This  can  be  understood  by 
considering  that  a  smooth  potential  variation  over  a  distance  of,  say,  1000,  would  require  a 

smooth  V2(|)  which  is,  in  turn,  given  everywhere  by  the  charge  density.  Maintaining  a  smooth 
charge  density  distribution  is  difficult  when  any  errors  in  determining  (ne-nO  are  multiplied  by  a 
huge  A'  .  There  is  one  effective  remedy  to  this  dilemma  described  in  Parker  ,  but  the  process 
reported  here  appears  to  be  more  efficient  in  the  short  Debye  length  limit.  This  method  involved 

TO 

the  combination  of  two  concepts.  One  uses  a  partial  implicitization  of  the  repelled  density  .  The 
other  simply  reduces  the  charge  density  to  an  acceptable  level  whenever  the  first  method  is 
inadequate.  An  “acceptable  level”  means  a  value  not  so  great  that  it  will  cause  spatial  oscillations 
in  the  potential. 

3.6.1  Implicitization 

The  right  hand  side  of  Equation  (158)  is  the  dimensionless  charge  density. 

q(r,®v(r))=A-2(ni(®v)-ne(®v))  (161) 

The  charge  density  at  the  present  iteration  may  be  linearized  about  the  previous  potential 
iteration 


q(®v)  =  q(®v  1)  +  q'(®vl)(®v -®v  ’) 


<3q 

where  q'  =  — - ,  and  the  r  dependence  has  been  dropped  for  clarity.  With  this  expression,  we 
0® 

may  write  the  implicit  Poisson  iteration  scheme 

-V2®v  -q'(®v~1)®v  =  q(®v~1)-q'(®v~1)®v_1  ^62 

Though  it  is  not  immediately  obvious,  the  implicit  character  of  Equation  (162)  makes  it  more 
stable  than  the  scheme  expressed  in  Equation  (160)  provided  q'  is  constrained  to  be  non¬ 
negative.  This  can  be  understood  by  realizing  that  in  Equation  (160)  the  charge  density  is  treated 
as  an  independent  variable,  whereas  in  Equation  (162)  the  charge  density  is  determined 
simultaneously  with  the  potential. 

The  finite  element  approximation  to  Equation  (162)  produces  the  matrix  equation 
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(163) 


=S-S'e®v_1 

e 

where  (e)  refers  to  each  element,  W  (which  produces  the  square  of  the  electric  field)  corresponds 
to  the  Laplacian  operator,  V  (which  produces  the  square  of  the  potential)  corresponds  to  the 
screening  matrix,  and  the  coefficient  of  the  screening  matrix,  S,  is  derived  from  q  by  the 
following  analysis. 

3.6.2  Charge  Limiting 

For  A>  1,  S  is  simply  the  total  charge  associated  with  each  node,  Q.  (Because  the  Debye  length 
is  longer  than  the  grid  resolution,  an  entire  grid  cube  of  plasma  with  one  species  eliminated  alters 
the  potential  by  no  more  than  the  temperature.)  However,  for  A  «  1,  numerical  noise  and 
features  like  a  sheath  edge,  which  may  span  only  a  few  XD,  become  incorrectly  amplified  when 
the  q  detennined  at  a  point  becomes  multiplied  by  the  entire  nodal  volume.  When  it  is  not 
possible  to  reduce  the  volume  element  size,  stability  can  be  preserved  by  replacing  Q  (and  Q ') 
with  a  reduced  value  S,  (and  S')  which  is  limited  by  the  maximum  allowable  charge  for  the 
element. 

Because  of  the  artificial  amplification  argument,  S  is  often  the  more  realistic  total  for  an  element, 
in  the  sense  that  it  produces  potential  variation  appropriate  to  the  spatial  resolution  rather  than 
causing  unphysical  overshoots.  Before  deriving  S,  we  define  the  barometric  potential  ®b  =  Inin,), 
which  (in  cases  where  electron  density  is  governed  by  the  Boltzmann  factor,  e°)  is  the  potential 
for  which  Q  =  0.  (For  the  “Linear”  and  “Nonlinear”  space  charge  density  fonnulations,  the  ion 
charge  density  is  also  dependent  on  the  potential.  For  these  cases,  the  barometric  potential  as 
used  here  is  zero,  ®b  =  0.)  Note  that  it  is  important  that  S— >  Q  as  ®  — >  ®b  if  quasineutral  regions 
are  to  be  modeled  correctly.  To  determine  S,  consider  a  capacitor  with  potential  difference 
®b  -  <f>,  area  L“,  and  a  separation  of  L.  The  charge  qc  on  this  capacitor  is  given  by 

p  1 2 

Qc  =  CAV  =  — ^ — (®b  -  ®)e0  (164) 

Le 

We  could  then  set  the  maximum  allowable  charge  per  element  to  be 

fiLimit  =a(®b-®)  (165) 

with  the  parameter  a,  set  to  the  maximum  value  consistent  with  the  stability  of  the  Poisson  solver. 
Thus  at  each  node,  the  charge  would  be 


S|  =  rninjqLimit|,|q|) 


(166) 


with 
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(167) 


-a, 


S'  = 


5q 


1 50’ 


for  S  =  q  Limit 


for  S  =  q 


However,  this  algorithm  can  easily  give  a  discontinuous  value  for  S which  could  lead  to 
numeric  instabilities. 


In  practice  a  different  formulation  is  used.  If  a  problem  has  been  specified  where  a  boundary 
potential  would  be  screened  in  less  than  an  element  or  two  (the  limit  of  any  code’s  resolution), 
sufficient  sheath  charge  is  redistributed  to  allow  the  potential  to  be  screened  over  the  minimum 
number  of  elements  that  is  consistent  with  stability.  Consider  the  solution  of  the  Helmholtz 
equation  in  a  single,  linear  zone.  The  variational  principle  for  the  Helmholtz  equation  is 


0  =  — | 

f 

f*  - 

1 

dV 

5® 

[J 

(V®)2  + 


(168) 


The  squared-gradient  term  is  minimized  by  a  constant  potential,  whereas  the  squared-potential 
tenn  is  minimized  by  going  symmetrically  through  zero.  So,  there  must  be  some  minimum  value 
of  A  for  which  the  solution  does  not  change  sign.  Minimizing  the  variational  function  on  the 
interval  [0,  1]  with  respect  to  ®  (1)  when  ®  (x)  =  (1-x)  ®  (0)+x  ®  (1)  gives 


0(1)  =  0(0) 


2- 


3A/ 


2  +  - 


3A" 


(169) 


9 

so  that  the  condition  for  non-oscillatory  potentials  becomes  6A  >  1 .  Thus  if  the  charge  is 
artificially  constrained  so  that  the  Debye  length  is  greater  than  V6  L  =  2.45 L  ,  the  Poisson 
solution  is  stable. 


For  each  charge  density  formulation,  a  computational  temperature  is  defined  and  used  so  that  the 

0D  8  d) 

Debye  length  analog,  A,  defined  by  — —  =  — ,  is  no  less  than  a  user  specified  fraction  (usually 

5(j)  A 

'A)  of  the  local  mesh  spacing.  The  formulas  for  S  and  S'  for  each  charge  density  fonnulation 
appear  in  Section  3.6.5. 


3.6.3  Analysis  of  the  Space  Charge  Stabilized  Poisson  Method 

The  charge  stabilized  Poisson  method  calculates  for  each  node  the  maximum  allowable  charge 
that  is  consistent  with  the  stability  of  a  linearly  interpolating  Poisson  solver.  This  method  is 
developed  above,  but  a  further  analysis  is  presented  here  to  help  the  user  interpret  its  impact. 

Nascap-2kA  charge  stabilization  is  accomplished  through  the  process  of  charge  limiting, 
illustrated  in  Figure  2 1 .  This  figure  shows  two  charge-versus-potential  curves  for  the  case  where 
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the  ion  density  is  fixed  (tracked  or  equal  to  the  background  density)  and  the  electron  density  is 
barometric.  Equation  (162)  is  rewritten  as 


q  =  aexp (-Om) (exp Ob -expO)  (1?0) 

where  Om  -ln(ak2Debye/L2)  a  is  user  specified  and  Ob  is  the  barometric  potential,  Ob  =  1  n (  n , ) .  For 

curve  qi,  Ob  =  -3.0  (n;  =  0.05);  for  curve  q2,  Ob  =  0.0  (n;  =  1.0)  and  for  both  curves  Om  =  -2.2, 
and  a  =  1.  For  each  curve,  the  limiting  charge  as  given  by  Equation  (165)  is  also  shown.  The 
limiting  charge  is  rewritten  here  as 

qumit  =  a  (Ob -O)  (171) 


which  intersects  the  “natural”  charge  curves  at  Oc  and  Ob.  The  charge  stabilization  method 
reduces  the  charge  to  the  limiting  value  when  O  >  Oc  and  uses  the  natural  charge  for  O  <  Oc. 

The  parameter  Om  provides  a  good  measure  of  the  limiting  process.  From  Figure  21  it  can  be 
seen  that  Om  is  the  point  at  which  the  slope  of  the  natural  charge  curve  equals  that  of  the  limiting 
charge  line.  Figure  22  shows  a  family  of  curves  giving  the  dependence  of  the  cutoff  potential  Oc 
on  the  barometric  potential  for  various  values  of  Om.  These  curves  were  obtained  by  numerically 
solving  for  the  zeros  of  the  difference  between  q  and  qLimit.  This  difference  equation  always  has 
two  solutions,  one  at  cpc  and  one  at  <T>b,  with  the  exception  of  a  degeneracy  at  cpc  =  cpb  =  <T>m, 
which  is  indicated  in  the  figure.  This  figure  shows  that  the  charge  limiting  is  minimal  for  cpm  >  -1, 
and  quite  severe  for  <f>m  <  -6  or  so. 


Figure  21.  Plots  of  Space  Charge  (Curves  qi  and  q2)  as  a  Function  of  Potential  as  Given  by 
Equation  (162).  The  Straight  Lines  Represent  the  Maximum  Allowable  Charge  for  Non-oscillatory 
Potentials.  The  “Natural”  Space  Charge,  qi  or  q2  is  Acceptable  for  Cases  for  Which  Slopes  of  the 
Curves  and  the  Corresponding  Line  Are  Equal. 
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Figure  22.  Plot  of  the  Space  Charge  Cutoff  Potential,  ®„  Versus  Barometric  Potential  (<Db  =  In  nd 
for  a  Series  of  ®m  Values  (-02,  -0.5,  -1.0,  -2.0,  -3.0,  --4.0  ...  -11.0).  The  Point  at  Which  ®m  =  <Db  =  ®c 
is  Also  Indicated. 

3.6.4  Sheath  Boundary  Potential 

Space  charge  limiting  must  be  accounted  for  when  choosing  the  sheath  boundary  potential. 
Space  charge  limiting  limits  the  rate  at  which  space  charge  can  cause  the  potential  to  drop. 

Consider  the  first  volume  element  of  a  sheath  to  satisfy  the  laws  of  Child  and  Langmuir  (planar 
space  charge  limiting).  At  the  sheath  edge  (z=0)  and  one  element  in  (z  =  L)  the  potential  and 
electric  field  are  given  in  Table  1 . 

Table  1.  Potential  and  Electric  Field  Variation  Given  by  Planar  Space  Charge  Limiting. 


Position 

L 

Potential 

K(L)4/3 

Electric  field 

(4/3)K(L)1/3 

By  Gauss’s  law,  the  charge  per  unit  area  in  this  element  is  given  by 

-L  =  1ic(lF 

£da  3  (i72) 
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2 

Nascap-2k  computes  the  charge  density  to  be  limited  by  (a/L~)  times  the  mean  potential 
(assuming  linear  interpolation),  and  so  gets 


Q 


s  A 


(173) 


Equating  these  two  expression  gives  a  =  8/3. 

Because  of  the  requirements  of  disk  space  and  computational  time  of  three-dimensional  codes, 
Nascap-2k  is  frequently  operated  at  high  Om  values,  i.e.,  a  coarse  mesh  with  respect  to  the  Debye 
length.  In  these  cases  charge  is  removed  at  almost  all  points.  Ideally,  the  charge  that  was 
removed  was  excess  charge  generated  by  the  coarse  gridding.  This  is  the  artificial  charge 
amplification  argument.  However,  since  Nascap-2k  must  be  reliably  stable,  the  result  is  that  too 
much  charge  is  often  removed.  This  results  in  an  enlarged  sheath  thickness  for  high  negative  Orn 
problems.  In  order  to  compensate  for  this  sheath  enlargement,  the  sheath  boundary  potential  must 
be  adjusted. 

Equating  the  above  equation  to  the  code  resolution,  L,  gives 

<(>x  =  5.1  x  lO“6L4/301/3n2/3  =  O.74(L/kD)4/30  n?zn 


The  potential  cj)x  may  be  interpreted  as  the  potential  below  which  Nascap-2k  underestimates 
screening.  At  best,  beyond  the  <j)x  contour,  the  potential  drops  about  one  order  of  magnitude  per 
element.  For  9  =  0. 1  eV,  n  =  10  m‘ ,  and  L  =  0.2  m,  we  find  §x  =  6  V.  If  the  6  V  contour  is  correctly 
placed,  the  0.6  V  contour  lies  at  least  one  element  beyond  (at  the  approximate  sheath  location),  and 
the  default  sheath  contour  (0.07  V)  is  yet  another  element  farther.  This  would  produce  a  sheath  area 
that  is  too  large.  The  suggested  criterion  for  the  sheath  boundary  potential  is 

^sb  =  Max(01n2,0.24c^x)  (175) 

where  0.24  =  e'L0/0'7  is  the  planar  screening  per  element  allowed  by  Nascap-2k.  Note,  however, 
that  <j)r  depends  strongly  on  the  grid  in  which  the  sheath  is  found,  so  that  if  an  increase  in  object 
potential  moves  the  sheath  from  grid  3  to  grid  2,  a  corresponding  larger  value  of  c()sb  should  be  used. 

3.6.5  Charge  Density  and  Derivative  in  Nascap-2k 

Nascap-2k  has  a  selection  of  eight  different  expressions  to  describe  the  charge  density  for 
different  kinds  of  problems.  Each  expression  uses  implicitization  and  charge  limiting  to  a 
different  extent.  As  there  is  no  space  charge  for  the  Laplace  space  charge  option,  the  following 
discussion  does  not  apply  to  that  case.  The  expressions  used  for  the  space  charge  and  its 
derivative  including  limiting  are  listed  below. 

The  parameter  D  is  the  local  mesh  spacing,  L,  divided  by  the  user  provided  Debye  limiting 
value,  which  defaults  to  2.  The  parameter  g  is  the  maximum  of  plasma  density  reduction  factor 
computed  by  neutral  wake  model  (0<g<  1 )  and  10‘6. 
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Laplace.  The  Laplacian  space  charge  option  specifies  that  the  charge  density  is  zero  and  no 
space  charge  limiting  is  needed. 


d(p/So)_Q 

d(j) 


i.e.,  charge  exists  only  on  object  surfaces  and  external  boundaries,  as  detennined  by  the 
boundary  conditions.  Space  charge  iterations  may  still  be  required,  however,  for  numerical 
reasons. 


(176) 


Linear  (Debye  Shielding).  The  Linear  space  charge  option  solves  the  Helmholtz  or  Debye  - 
Huckel  equation: 


P_  =  _A_ 

£0  Ki 

Xnl2  =  max(A,p  /g,D2) 

d(p/s0)  =  1 

#  “x2nl 


(177) 


Nonlinear.  Nonlinear  space  charge  is  appropriate  for  most  low-Earth-orbit  type  plasmas  in  the 
presence  of  high  applied  potentials.  Poisson’s  equation  is  solved  with  space  charge  given  by: 


p/s0 


max(l,C((j),E)) 


C  (<)>,  E )  =  min  (( Rsll /r  )2 , 3 . 545  |4,/enl  L 2 ) 
(R„,/r)2  =2.29|ex„./©„]|,2“  |e„,/<t>|“509 

K\=  max(A,p  /g,D2) 

0m  =  eftig  /K)m 


(178) 


The  term  C(<|>E)  (analytic  focusing)  comes  from  fitting  a  finite  temperature  spherical  (Langmuir- 
Blodgett)  sheath.  If  analytic  convergence  is  turned  off  (on  the  Potentials  Advanced  screen)  then 
C(c()  E)  is  set  to  zero.  Additional  adjustments  are  made  to  the  convergence  portion  of  this  formula 
for  significant  values  of  velocity. 
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Note  that  the  derivative  of  the  convergence  contribution  to  the  charge  density  with  respect  to  the 
potential  is  assumed  to  be  negligible. 

Frozen  Ion.  The  Frozen  Ion  formulation  is  intended  for  short  timescale  (typically  sub¬ 
microsecond)  problems  for  which  it  is  a  good  approximation  to  assume  that  ions  remain 
stationary  and  at  ambient  density  (“ion  matrix”  approximation),  but  electrons  achieve  barometric 
equilibrium.  The  space  charge  function  depends  on  the  mesh-dependent  potential,  (jq  <  0,  which 
satisfies 


i-exp(cj)1/e)  =  -(^/D)2(())1/e)  (180) 

for  D>a.  Otherwise,  (j)i=-lxl0"6.  The  meaning  of  this  quantity  is  that  the  charge  density  increases 
linearly  at  a  stable  rate  for  negative  potentials  between  zero  and  (jq,  and  then  asymptotes  to  the 
ambient  ion  density  for  potentials  more  negative  than  (jq. 

The  space  charge  is  then  given  by 
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0  >  <j)  >  <j)j 


(181) 


(182) 


The  derivative  of  charge  with  respect  to  potential  appears  in  the  variational  function  as  a 
coefficient  of  cj>  .  If  this  is  too  large,  then  minimizing  (j)  means  you  oscillate  from  one  grid  point 
to  the  next  in  order  to  get  lots  of  intennediate  zeroes,  which  is  not  a  satisfactory  solution. 

Barometric  .  The  Barometric  algorithm  is  appropriate  for  cases  in  which  all  the  surfaces  are  at 
potentials  comparable  to  or  below  the  plasma  temperature  (or  negative)  and  there  is  a  region  of 
low  density,  such  as  a  plasma  wake.  The  ion  density  is  taken  as  the  plasma  density  decreased  by 
the  wake  factor  g.  If  the  ions  are  sufficiently  dense,  quasi-neutrality  requires  the  potential  to  be 
the  “barometric”  potential,  (pb  (defined  below).  However,  because  the  plasma  may  be  under- 
dense  or  influenced  by  nearby  surfaces,  the  electron  density  must  be  allowed  to  differ  from  the 
ion  density.  We  take  the  electron  density  to  be  given  by  the  expression 

P^  =  _  p|^  fexp(min(1°,((t>-(t>b)/eb))  for  (j)  -  (j)b  <  0 
so  80  j  l  +  min(l0,(<|>-<|>b)/eb)  for  <|>  -  <|>b  >  0 
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Full  Trajectory  Ions.  Ion  densities  are  calculated  from  steady-state  ion  trajectories.  Electrons 
are  barometric. 
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(186) 


Hybrid  PIC.  This  algorithm  is  used  for  timescales  (typically  sub-millisecond)  for  which  it  is 
practical  to  treat  ion  motion,  but  electrons  may  be  considered  in  barometric  equilibrium.  The  ion 
density  is  computed  from  actual  ion  macroparticles.  The  electron  charge  density  is  barometric  for 
negative  potentials  and  increases  linearly  (orbit  limited)  for  positive  potentials. 

Pe___pJ_J  exp(min(lO,((|)-4)by0hP)) 

e0  g0  exp  (min(l  0  ,-4>b  /  0bp  )j(l  +  4>/0hp ) 


for  (j)  <  0 
for  0  <  4> 


(187) 
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When  the  potential  is  near  the  barometric  potential  (|cj)  —  (j)b  |  <0!),  space  charge  limiting  is  not 

needed,  and  0hp  is  0.  When  the  potential  is  far  from  the  barometric  potential  (02  <  |c()-c()b|),  the 

effective  temperature  used  is  the  space  charge  limited  value.  At  intermediate  potentials,  the 
effective  temperature  varies  linearly. 
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hip 


is  the  value  of  0hP  at  (j)  =  0  and  is  only  needed  for  positive  potentials. 
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Full  PIC.  For  this  option,  it  is  assumed  that  the  total  charged  density  (electron  and  ion)  was 
previously  computed  and  stored. 
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3.7  Macroparticles 

Macroparticles  represent  either  current  or  charge.  Macroparticles  that  represent  current  are  used 
for  steady-state  problems  and  ones  that  represent  charge  are  used  for  dynamic  problems. 
Macroparticles  are  also  used  for  visualization. 

3.7.1  Macroparticle  Creation 

When  macroparticles  are  created  at  a  sheath  surface,  the  sheath  surface  is  divided  into  triangles 
and  a  macroparticle  representing  the  current  from  that  section  of  the  sheath  is  created  at  the 
center  of  the  triangle. 

When  macroparticles  are  created  at  the  problem  boundary,  the  outer  surface  of  each  grid 
boundary  volume  element  is  divided  into  squares  (of  size  requested  by  the  user)  and  a 
macroparticle  is  created  at  the  center  of  each  square  to  represent  either  the  current  through  that 
area  (“Boundary”  and  “BField”)  or  the  charge  passing  through  that  area  during  a  user  specified 
time  (“Boundary  Injection”). 

When  macroparticles  are  created  throughout  the  volume  (for  uniform  charge  density 
initialization  or  charge  exchange),  they  are  created  at  points  appropriate  to  the  weighting 
functions. 

When  macroparticles  are  created  at  user  specified  locations  with  user  specified  velocities,  each 
macroparticle  represents  current. 
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When  macroparticles  are  created  on  surface  elements,  the  surface  is  divided  into  the  requested 
number  of  quadrilaterals  and  triangles  and  macroparticles  are  created  at  the  center  of  each.  The 
macroparticles  are  distributed  evenly  in  velocity,  azimuthal  angle,  and  cosine  of  the  polar  angle. 

For  visualization  only,  macroparticles  can  also  be  created  at  the  intersection  of  the  sheath  surface 
and  a  user  specified  plane. 

3.7.2  Macroparticle  Tracking 

Macroparticles  are  tracked  using  a  third  order  energy  conserving  algorithm. 

x(t  +  dt)  =  x(t)+v(t)dt  +  -^-dt2  +  — dt3  +  — dt4 
v  ’  v  ’  v  7  2  6  24 

v(t  +  dt)=  v(t)+x2dt  +  — dt2  +  — dt3 

2  6  (192) 


where 


x2  =—  (E  +  vxB) 
m 

x3  =  —  (v  •  VE  +  x,xB) 
m 

x4  =-(x,xB) 
m 


(193) 


The  electric  field  is  that  at  the  location  of  the  macroparticle  detennined  by  interpolation  from  the 
nodes  of  the  volume  element.  The  magnetic  field  is  the  sum  of  a  constant,  uniform  field  supplied 
by  the  user  and  the  sum  of  all  the  magnetic  dipoles  defined  during  object  definition.  As  the  field 
sources  are  not,  in  reality,  point  dipoles,  if  a  macroparticle  is  within  1  cm  (0.01  m)  of  a  dipole, 
that  dipole’s  contribution  to  the  particle’s  magnetic  field  is  calculated  as  if  it  were  1  cm  away.  If 
the  macroparticle  is  within  1  mm  (0.001  m)  of  a  dipole,  the  contribution  of  that  dipole  is  set  to 
zero. 


B(rl-B  +  V  3»(P  »)~P 

“V1/  “ambient  2-j  ,3 

all  dipoles  |r~ro|  (194) 

Here,  p  and  rQ  are  the  dipole  moment  and  location  and  n  is  a  unit  vector  from  r  to  r0, 

Macroparticles  are  tracked  for  a  timestep,  which  can  be  as  long  as  a  minute  for  static  problems  or 
as  short  as  microseconds  for  dynamic  problems.  The  timestep  is  divided  into  substeps  such  that 
the  macroparticle  can  travel  no  more  than  a  fraction  of  the  local  mesh  size  in  a  single  substep.  By 
default  the  fraction  is  0.1;  the  user  may  specify  an  alternative  value.  The  particle’s  total  energy  is 
calculated  at  the  beginning  of  the  timestep  and  preserved  during  the  timestep.  Changes  in 
potential  between  timesteps  are  reflected  by  a  change  in  the  particle’s  total  energy. 
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If  the  macroparticle  represents  charge,  at  the  end  of  each  timestep  and  optionally  at  the  end  of 
each  substep,  the  charge  and  current  are  assigned  to  the  nodes  of  the  element  the  macroparticle  is 
in,  using  the  same  weights  as  the  potential.  The  charge  density  (charge  divided  by  element 
volume)  is  assigned  to  the  element.  If  the  assignment  is  perfonned  at  the  end  of  each  substep,  the 
charge  assigned  is  divided  by  the  fraction  of  the  timestep  that  the  substep  represents. 

If  the  macroparticle  represent  currents,  at  the  end  of  each  substep  its  charge  density  (current 
times  substep  time)  is  optionally  assigned  to  the  element  within  whose  bounds  the  macroparticle 
is  located. 

When  a  macroparticle  impacts  a  surface  element,  its  current  or  charge  is  assigned  to  the  surface 
element. 

3.7.3  Macroparticle  Splitting 

Macroparticles  are  split  both  to  avoid  having  heavy  macroparticles  in  well-resolved  regions  and 
to  simulate  a  thennal  distribution.  Macroparticle  splitting  has  been  implemented  for 
macroparticles  read  from  an  external  file,  for  a  unifonn  distribution  of  macroparticles  for 
initialization,  and  for  macroparticles  injected  from  the  boundary.  Macroparticles  can  also  be  split 
when  entering  a  more  finely  resolved  grid. 

The  following  general  principles  are  used  in  generating  the  smaller  replacement  macroparticles. 

1.  Macroparticles  are  split  in  velocity  space  only.  Because  particles  frequently  need  to  be  split 
in  high-field  regions,  spatial  splitting  would  raise  problems  with  energy  conservation. 

2.  To  be  split  in  velocity  space,  a  macroparticle  must  carry  a  temperature.  We  assume  the 
temperature  is  always  isotropic.  The  fission  products  carry  half  the  temperature  of  the 
original  macroparticle,  while  the  remaining  thennal  energy  appears  as  kinetic  energy  of  the 
split  macroparticles. 

3.  For  splitting  purposes,  we  define  the  Z-axis  to  be  along  the  direction  of  the  macroparticle 
velocity,  the  X-axis  randomly  chosen  in  the  plane  normal  to  Z,  and  the  Y -axis  mutually 
perpendicular. 

4.  We  split  into  two  or  three  macroparticles  with  added  velocity  along  each  axis,  except  that  we 
may  elect  not  to  split  along  the  Z-direction  if  the  kinetic  energy  exceeds  the  thennal  energy. 
Not  splitting  along  Z  helps  ameliorate  macroparticle  proliferation,  but  makes  an  error  by  not 
preserving  the  original  macroparticle  temperature  along  Z.  We  thus  end  up  with  eight,  nine, 
or  twenty-seven  new  macroparticles. 

5.  Macroparticle  velocity  is  assumed  to  be  acquired  by  acceleration  rather  than  actual  drift  (i.e., 
spacecraft  velocity).  If  there  is  actual  drift  (e.g.,  ram  velocity),  then  the  drift  velocity  is 
removed  before  splitting  the  macroparticle,  and  added  back  after. 
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6.  If  the  macroparticle  with  speed  uo  is  split  by  two  along  the  X  or  Y  axis,  the  new  velocity  is 
±  0.707^/eG/m  .  Along  the  Z  axis,  the  velocity  increment  is  calculated  as  if  the  temperature 

f  i - 7 \ 


were  9-2mu; 


1  +  - 


0 


mu; 


■1 


7.  If  the  macroparticle  with  speed  uo  is  split  by  three  along  the  X  or  Y  axis,  there  is  a  zero- 
velocity  central  macroparticle  and  two  “probe”  macroparticles  with  velocity  ±  O.866^/e0/m 
Along  the  Z  axis,  the  velocity  increment  is  calculated  as  if  the  temperature  were 


0-2mu; 
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mu; 


-1 


Examples  of  the  effects  of  macroparticle  splitting  appear  in  Reference  39. 


3.7.4  Reverse  trajectory  technique  to  compute  currents  to  detectors 


In  the  reverse  trajectory  approach,  the  current  is  detennined  by  an  integral  over  the  thermal 
distribution  of  the  incident  charged  particles.  The  flux  to  a  specific  location  between  9min  and 
0max  of  the  normal  and  with  speeds  between  vrain  and  vmax  is  given  by  the  integral40 
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(195) 


where  v,  the  velocity  incident  on  the  detector,  is  expressed  in  spherical  coordinates  (v,  ,  tp) 

with  respect  to  the  surface  nonnal,  Voo  is  the  velocity  at  infinity,  and  U  is  the  ram  velocity. 
Portions  of  phase  space  for  which  Voo  does  not  exist,  do  not  contribute  to  the  integral.  The  Va> 
value  corresponding  to  a  v  exists  only  when  the  trajectory  of  the  particle  with  incident  velocity  v 
connects  with  infinity,  i.e.,  exits  the  computational  space.  Energy  conservation  relates  the 
magnitude  of  Voo  to  that  of  v, 


m  , 

— v  +  e4> 


(196) 


where  cj)  is  the  detector  surface  potential.  The  exponent  includes  a  dot  product,  and  is  therefore  a 
function  of  the  angle  between  the  velocity  vector  at  infinity  and  the  ram  vector.  It  is  a 
complicated  function  of  v  and  in  this  technique  is  detennined  by  tracking  macroparticles. 


When  written  as  a  sum,  the  above  integral  is 


f 

n 

v 


3/2 


Jm.x  coda^in  )  2u 

'y 1  y '  ^^(A(pAcos9Av)v2( 

v=vmillcos9=cos(9mal)  <p=0 


VCOS 
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This  quantity  is  computed  by  generating  and  tracking  a  distribution  of  macroparticles  with  initial 
velocities  distributed  in  each  of  v,  cos0,  and  (j>  and  with  weight 


n 


m 

27ie0 . 


.3/2 


(A(pAcos  0Av)v°  (vcos  If) .  The  contribution  of  each  macroparticle  is  the  product  of 


its  weight  and  exp 


m 


(vx  -  U)  |  where  v.y,  is  the  velocity  when  it  exits  the  computation 
v  2e0  J 


space.  Those  macroparticles  that  intersect  the  object  do  not  contribute. 

3.8  Particle-in-Cell  Tools 


3.8.1  Transverse  Surface  Currents 

We  developed  a  pseudopotential  approach  (similar  to  the  velocity  potential  used  to  describe 
potential  flow41  in  fluid  dynamics)  to  the  solution  of  the  current  continuity  equation  for  the 
computation  of  transverse  surface  currents  on  a  spacecraft  acting  as  an  antenna.  As  a  boundary 
condition,  one  surface  element  of  each  antenna  element  is  specified  as  connected  to  the  biasing 
power  supply,  and  the  solution  provides  the  required  current.  The  vector  transverse  surface 
current  in  each  surface  element  is  that  which  provides  the  best  fit  to  the  edge  currents. 

A  solution  having  no  circulating  currents  is  guaranteed  if  we  assign  a  pseudopotential  value  to 
each  element  of  surface  and  take  the  current  across  their  common  edge  to  be  proportional  to  the 
difference  in  their  pseudopotentials.  (In  that  case,  current  that  flows  “downhill”  would  need  to 
flow  back  “uphill”  in  order  to  circulate.)  The  current  is  also  taken  as  proportional  to  the  edge 
length  and  inversely  proportional  to  the  distance  between  the  centers  of  the  elements. 

The  result  of  the  above  treatment  is  a  matrix  that  multiplies  the  vector  of  pseudopotentials  to 
describe  the  buildup  of  charge  in  the  surface  elements,  which  is  then  set  equal  to  the  source 
vector  generated  from  Nascap-2k  results  at  each  timestep.  This  system  of  equations  is  then 
solved  using  the  ICCG  (Incomplete  Cholesky  Conjugate  Gradient)42  algorithm. 

Figure  23  shows  graphically  the  quantities  that  make  up  the  surface  current  equation  for  each 
surface  element.  Moving  the  surface  currents  to  the  left  side  of  the  equation,  and  all  other 
(known)  quantities  to  the  right,  and  expressing  the  edge  currents  as  proportional  to 
pseudopotential  differences,  results  in  a  sparse  symmetric  matrix  that  is  amenable  to  solution 
using  the  ICCG  algorithm  once  the  surface  element  connected  to  the  power  supply  (denoted  the 
“injection  element”)  is  set  to  a  fixed  pseudopotential.  The  current  required  to  bias  the  conductor 
may  be  obtained  by  evaluating  the  equation  associated  with  the  injection  element  and  can  be 
verified  to  be  the  sum  of  the  total  change  in  charge  less  plasma  currents  for  all  the  remaining 
surfaces  of  the  conductor.  Surface  currents  are  solved  separately  for  each  electrically  isolated 
component. 
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Total  change  in  charge  Plasma  Currents  Surface  Currents 

Figure  23.  The  Change  in  Charge  on  a  Surface  Element  is  Comprised  of  Plasma  Currents  and 
Surface  Currents. 


3.8.2  Propagating  Fields 

The  transverse  surface  currents,  the  volume  ion  currents  computed  during  particle  tracking,  and 
the  volume  electron  currents  computed  from  tracking  or  saved  in  the  database  by  an  external 
code,  are  a  source  of  propagating  electromagnetic  fields.  Nascap-2k  can  compute  the  magnetic 
field,  the  vector  potential,  and  the  rate  of  change  of  the  vector  potential  from  these  currents. 


A(r)=£ 
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P(r,  t)  =  s0  c2E(r,  t)x  B(r,  t) 


(199) 

(200) 

(201) 


where  the  sums  are  over  the  surface  elements  (i)  and  volume  elements  (k)  and  tm  is  the  time  at 
the  present  timestep,  and  P(r,t)  is  the  Poynting  vector. 

3.8.3  Electron  Currents  (Using  Pseudopotential) 

The  external  code  included  with  Nascap-2k  provides  a  non-PIC  (Particle-In-Cell)  method  to 
calculate  volume  electron  currents  that  satisfy  the  equation  of  continuity  and  are  consistent  with 
an  ion  dynamics  simulation.  The  volume  electron  current,  jeiec  is  taken  to  be  proportional  to  the 
gradient  of  a  “pseudopotential,”  vp:  jelec  =  oV\|/ ,  where  the  conductivity  tensor,  o,  depends  on  the 
electron  density  and  magnetic  field.  The  equation 


-  V  -<tV\|/ 


dP  elec 

dt 


(202) 


is  solved  by  minimizing  the  function 
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(203) 


f  (l 

dV  —  Vi|/  -oVvj/ -\\/ 
V2 


dpelec^ 

dt  , 


Like  the  solution  to  Poisson’s  equation,  the  integrals  are  calculated  for  each  volume  element.  The 
pseudopotential  is  taken  to  be  trilinear  within  each  volume  element,  so  it  is  given  by 


¥e(x)  =  XN.(xVi 


(204) 


where  “i”  are  the  nodes  of  the  element  indexed  by  “e”,  and  the  Ni  are  the  interpolants.  Within 
each  volume  element,  the  integral  is  given  by 


1  vpf  (n  \dNi(x)  dNj(x)^|  V  ff-\T  l  t^Pelec 

-S'fi'fj  X  aap(B>peiec)-r^--r —  dv-Xv i  Ni(x)— ff- 

21T  Jojl  dx«  dxP  J  i  jy  dt 

V,  V- 


dV  (205) 


As  a  and  peiec  are  taken  to  be  constant  within  each  volume  element,  the  integrals  depend  only  on 
the  volume  size  and  shape.  To  minimize  the  volume  integral  in  Equation  (205),  we  set  its 
derivative  with  respect  to  vp  equal  to  zero.  Thus  the  equation  to  be  solved  is 


XA^  =x 

ej  e 


dP  elec 

dt 


Nj(x)dV  , 

J 

Ve 


(206) 


with  the  finite  element  matrix  for  each  element  given  by 


^  afi 


dNi  dNj 
dx„  dx„ 


dV 


(207) 


The  576  (8  x  8  x  3  x  3)  matrix  elements  of 


PN,  5Nj 
0Xa  0Xp 


dV  are  precomputed  for  the  unit  cube.  For 


each  element,  Ay  is  an  8x8  matrix  depending  on  B  and  peiec-  Because  the  conductivity  tensor  is 
anisotropic  for  B40,  matrix  elements  that  are  zero  in  the  scalar  case  become  nonzero  in  the 
presence  of  magnetic  field. 


The  right  hand  side  is  formed  by  adding  one-eighth  of  the  charge  density  derivative  (formed  by 
differencing  the  present  electron  charge  density  with  that  at  the  previous  timestep)  for  each 
volume  element  to  each  of  its  nodes  after  multiplying  by  the  cube  of  the  mesh  size.  The  matrix 
elements  are  accumulated  using  a  sparse  matrix  storage  scheme  for  eventual  solution  using 


ICCG. 
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3.8.3. 1  Boundary  Conditions 

Inner  Boundary  Conditions  (Object  Currents) 

The  solution  is  not  extended  to  the  special  elements  (those  containing  surfaces).  Rather,  the 
following  boundary  condition  is  applied:  Electron  current  flows  through  the  boundary  between 
an  empty  and  a  special  cell  if  the  electric  field  points  away  from  the  special  cell,  and  no  electron 
current  flows  in  the  opposite  polarity. 

If  an  empty  volume  element  neighbors  a  special  element  and  has  its  electric  field  pointing  away 
from  the  special  element,  then  electrons  are  assumed  to  flow  across  the  interface  to  the  special 
element  at  the  local  plasma  thermal  current.  Each  of  the  four  nodes  on  the  corresponding  face 

2  P 

jth  ^£- ,  where  j*  is  the  electron  plasma  current  in 
Pa 

the  ambient  plasma  and  pa  is  the  ambient  plasma  density.  If  the  electric  field  points  toward  the 
special  element,  then  no  electron  current  crosses  the  interface,  and  no  action  is  taken. 


has  its  right  hand  side  augmented  by  -^-(Ax) 


Outer  Boundary  Conditions 


The  pseudopotential  is  assumed  to  vary  inversely  with  radius  (measured  from  the  grid  center) 
outside  the  computational  space.  As  a  consequence,  for  each  outer  boundary  square,  each  of  its 

nodes  is  assigned  an  additional  diagonal  matrix  element  given  by  ^  —  oaprarp  ,  where  Q  is 

ap 

the  solid  angle  subtended  by  the  outer  boundary  square  relative  to  the  grid  center,  r  is  the  vector 
from  the  grid  center,  and  a  is  the  conductivity  tensor. 


3. 8.3.2  Grid  Interface 

Grid  interfaces  are  constrained  to  have  the  resolution  of  the  coarser  grid.  The  dislocations  at  grid 
interfaces  can  make  convergence  to  a  stable  solution  difficult. 

Since  the  matrix  elements  for  the  entire  computational  grid  are  accumulated  in  order  to  use  an 
ICCG  solver,  the  matrix  must  be  modified  to  place  constraints  on  inner  grid  boundary  nodes  and 
assign  matrix  elements  generated  in  the  boundary  volumes  of  the  inner  grid  to  outer  grid  nodes. 


3. 8.3.3  Upwind  Correction 

An  upwind  technique  is  used  to  reduce  discontinuities  at  grid  interfaces.  The  input  current  to 
each  element  is  calculated,  and  then  the  magnitude  of  the  cell  current  (now  redefined  as  output 
current)  is  adjusted  (after  accounting  for  the  known  charge  accumulation  rate)  to  maintain  charge 
conservation.  The  current  direction  calculated  by  the  finite  element  solution  is  preserved.  The 
following  algorithm  is  used. 

Conservation  of  charge  in  volume  element  i  can  be  written  as  (lin  -  Iout  )At  =  Aq ,  where  Im  is  the 

current  of  electrons  flowing  into  element  i  from  neighboring  volume  elements,  Iout  is  the  current 
of  electrons  flowing  from  element  i  to  neighboring  volume  elements,  and  Aq  is  the  change  in 
electron  charge  observed  in  the  simulation:  Iin  =  ^  Jj  'nii^ij  >  where  J,  is  the  calculated 

i3  J  j  'nij  >o 


75 

Approved  for  public  release;  distribution  is  unlimited. 


electron  current  density  in  neighboring  cell  j,  ny  is  the  outward  normal  from  cell  i  into  cell  j,  Ay 
is  the  interface  area  between  the  two  cells,  and  the  inequality  restricts  the  sum  to  volume 
elements  j  from  which  electrons  flow  into  element  i.  Similarly,  Iout  =  -  ^  J,  •  n  y  A  y  ,  where 

jaJj-n^O 


the  inequality  now  restricts  the  sum  to  volume  elements  j  to  which  electrons  flow  from  element  i. 
We  can  make  conservation  of  charge  true  by  scaling  the  current  in  volume  element  i: 


J i — ^  Ji 


linAt- Aq 


I 


out 


.  However,  this  affects  the  conservation  of  charge  in  “downwind”  cells,  so  an 


iterative  “relaxation”  process  must  be  implemented  to  achieve  an  improved  solution. 


Five  iterations  of  this  algorithm  are  applied. 
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LIST  OF  SYMBOLS 


The  symbols  in  the  following  table  are  used  in  this  document. 


bj 

Range  parameters 

d 

Thickness  or  distance 

e 

Magnitude  of  the  electron  charge  (1.60  x  10'19  C) 

f 

Distribution  function 

g 

Reduction  in  plasma  density  in  wake  in  the  absence  of  surface  potentials 

i  and  j  and 
k  and  m 

Indices 

Jth 

Plasma  thermal  current 

jx 

Component  of  the  current 

l 

Length 

m 

Particle  mass 

n 

Component  density 

na 

Ambient  plasma  density 

n 

Surface  normal  or  unit  vector 

p 

Magnetic  moment 

q 

Charge  on  a  surface  or  particle  charge 

q; 

Range  parameters 

r 

Distance  or  location  in  space 

ro 

Location  of  magnetic  moment 

s 

Axes  of  sun  pointing  coordinate  system 

t 

Time 

v  or  v 

Velocity 

u  or  v 

Parameters  in  BEM 

x  or  x 

Position  or  path  length 

y 

Substitute  integrand 
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A 

Area  or  Point  of  triangle  in  BEM  discussion 

A 

BEM  or  pseudopotential  matrix 

B  and  B 

Magnetic  field  and  magnetic  field  magnitude  or  point  of  triangle  in  BEM  discussion 

B  max 

Estimate  of  maximum  possible  barrier  height  to  the  escape  of  photoelectrons 

c 

Capacitance  or  point  of  triangle  in  BEM  discussion  or  Convergence  factor  or  coefficient 

DCl 

Child  Langmuir  distance 

D 

Electric  field  (in  Maxwell  equation) 

E 

Particle  energy  in  eV 

E  and  E 

Electric  field  and  electric  field  magnitude 

E0 

Parameter  of  Gaussian  component  of  Fontheim  distribution  function 

Ei 

Barrier  height  to  escape  of  photoelectrons 

F 

■^max 

Energy  at  which  secondary  yield  due  to  normally  incident  electrons  peaks 

F 

Flux  or  interpolant 

F 

BEM  matrix 

G 

BEM  matrix 

H 

Heaviside  step  function 

G 

Interpolant 

I 

Current  (amperes) 

I 

Current  vector 

Isun 

Relative  sun  intensity 

j 

Current  density  vector 

L 

Mesh  spacing  or  lower  limit  of  energy  or  velocity  integral 

N 

Number  of  nodes  of  surface  element 

Ni 

Interpolant 

P  or  P 

Point  of  triangle  in  BEM  discussion  or  point  in  volume  or  Poynting  vector 

Q 

Finite  element  matrix 
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R 

Resistance  or  Sphere  radius  or  Escape  depth  (range)  or  point  in  volume 

S 

Surface  area  or  Stopping  power  or  Surface  index  or  reduced  charge  density 

AT 

Time  interval 

T 

Pairwise  scalar  products  in  BEM 

U 

Velocity  of  plasma  with  respect  to  spacecraft,  -Vsc 

V 

Volume  or  BEM  coefficients 

vsc 

Spacecraft  velocity 

w 

Finite  element  matrix 

Y 

Yield 

z 

Atomic  number 

a 

Parameter  used  in  ion  secondary  yield  or  Parameters  used  in  space  charge  limiting  or  parameter  used 
to  specify  power  law  component  of  Fontheim  distribution  function 

p 

Parameter  used  in  ion  secondary  yield 

X 

Angle  between  flow  vector  and  the  particle  velocity 

^max 

Peak  yield  for  normal  incidence 

Permittivity  of  vacuum  (8.854  x  10'12  F  m"1) 

8 

Relative  dielectric  constant 

<l> 

Surface  potential  or  Space  potential 

Barometric  potential 

Backscatter  coefficient  or  Fraction  of  the  distribution 

K 

Parameter  used  in  kappa  distribution  function 

XD 

Debye  length 

V 

Iteration  index 

e 

Plasma  temperature  in  eV 

p 

Resistivity  or  Charge  density 

a  or  a 

Surface  charge  density  or  conductivity  tensor 

82 

Approved  for  public  release;  distribution  is  unlimited. 


C?CEX 

Charge  exchange  cross  section 

5 

Small  number  used  to  evaluate  current  integral 

V 

Incident  angle  or  Angle  between  neutralizer  axis  and  look  direction  or  pseudopotential 

c 

Coefficient  in  Fontheim  formula 

^gauss 

Parameter  of  Gaussian  component  of  Fontheim  distribution  function 

A(j) 

Difference  in  potential  across  dielectric  layer 

0> 

Potential  vector 

Q 

Solid  angle 
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